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Abstract 

Functionally graded materials (FGMs) have become one of the major research topics in the 
mechanics of materials community during the past fifteen years. FGMs are heterogeneous ma- 
terials, characterized by spatially variable microstructures, and thus spatially variable macro- 
scopic properties, introduced to enhance material or structural performance. The spatially 
variable material properties make FGMs challenging to analyze. The review of the various 
techniques employed to analyze the thermomechanical response of FGMs reveals two distinct 
and fundamentally different computational strategies, called uncoupled macromechanical and 
coupled micro-macromechanical approaches by some investigators. The uncoupled macrome- 
chanical approaches ignore the effect of microstructural gradation by employing specific spatial 
variations of material properties, which are either assumed or obtained by local homogenization, 
thereby resulting in erroneous results under certain circumstances. In contrast, the coupled ap- 
proaches explicitly account for the micro-macrostructural interaction, albeit at a significantly 
higher computational cost. The higher-order theory for functionally graded materials (HOT- 
FGM) developed by Aboudi et al. (1999) is representative of the coupled approach. However, 
despite its demonstrated utility in applications where micro-macrostructural coupling effects are 
important, the theory’s full potential is yet to be realized because the original formulation of 
HOTFGM is computationally intensive. This, in turn, limits the size of problems that can be 
solved due to the large number of equations required to mimic realistic material microstruc- 
tures. Therefore, a basis for an efficient reformulation of HOTFGM, referred to as user-friendly 
formulation, is developed herein, and subsequently employed in the construction of the efficient 
reformulation using the local/global conductivity matrix approach. In order to extend HOT- 
FGM’s range of applicability, spatially variable thermal conductivity capability at the local level 
is incorporated into the efficient reformulation. Analytical solutions to validate both the user- 
friendly and efficient reformulations are also developed. Volume discretization sensitivity and 
validation studies, as well as a practical application of the developed efficient reformulation, are 
subsequently carried out. The presented results illustrate the accuracy and implementability of 
both the user-friendly formulation and the efficient reformulation of HOTFGM . 


1 Introduction 

Functionally graded materials (FGMs) have become one of the major research topics in the me- 
chanics and materials communities during the past fifteen years. FGMs are heterogeneous materials 
which are composed of two or more phases with different properties, and therefore they belong to 
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the broad class of materials known as composites. Unlike traditional composites, however, the spa- 
tial distribution of phases in FGMs is not statistically homogeneous. As the name suggests, FGMs 
are materials with macroscopic properties characterized by spatial gradients. This is accomplished 
by using variable spacings between individual inclusions or by using inclusions with different prop- 
erties, sizes and shapes. Research results obtained thus far have demonstrated that FGMs have 
great potential for improving material/structural performance in many engineering applications 
precisely because of their spatially graded heterogeneous microstructure. In particular, grading or 
tailoring the internal microstructure of a composite material or structural component enables the 
designer to integrate both the material and structural aspects into the final design and product. 
This offers a number of advantages over traditional methods of tailoring composite materials and 
structural components through the elimination of sharp material discontinuities which often lead 
to undesirable stress concentrations, thereby degrading the structural component’s strength and 
durability. 

The concept of modern man-made FGMs was first proposed in 1984 by material scientists in 
the Sendai area of Japan as a means of developing thermal protection materials. In particular, this 
concept was originally conceived for thermal protection of fuselage exterior and engine materials 
intended for use in high-speed civil transport aircraft under development in the mid 1980’s through 
mid 1990’s, and then rapidly spread to other fields including physics, chemistry, agriculture, biology 
and medicine. After almost twenty years of development, the term FGM is now widely used in 
various fields to describe this class of materials. Strictly speaking, however, graded materials are not 
new. They can be found in nature, and they also have been utilized extensively by human beings 
for many years without being explicitly called FGMs. The examples include graded materials 
developed long ago, such as case-hardened steel, which are still in common use today, Kaysser and 
Ilschner (1995). Other naturally occurring biological FGMs include bamboo, which has been used 
for a long time for structural components, corn, barley, etc., Figure 1. Nevertheless, the new and 
exciting thing about modern FGMs is the realization that the property gradient can be designed 
at the microstructural level to tailor the material’s specific function and performance required by 
engineering needs. 





(a) Bamboo | 1Qmm | (b) Corn Q-2 mnn ( C ) Barley 0.2 mm 

Cross Section 

Figure 1. Selected FGMs found in nature. (Source: Amada, 1995). 


As stated in the foregoing, FGMs are excellent candidates for applications involving severe 
thermal gradient environments. Some practical examples of the successful application of the FGM 
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concept include zirconia coatings deposited on a nickel-chromium alloy substrate by plasma spray- 
ing, where the composition profiles of the graded interlayers are carefully optimized to minimize 
thermally-induced stresses under service loading conditions. By appropriately adjusting the mi- 
crostructural transition, that is by placing the ceramic-rich part of a functionally graded material 
in the high temperature environment and the metallic-rich part in the low temperature environ- 
ment, with a gradual microstructural transition in the direction of the temperature gradient, it 
is possible to achieve optimum temperature, deformation and stress distributions. An example 
of a graded thermal barrier coating is given in Figure 2. In this case, the requirements for heat- 
resistant and oxidation-resistant properties are fulfilled, as well as those for adequate mechanical 
toughness. In another successful application, Honda Engineering Co. Ltd. developed a functionally 
graded cemented carbide with improved wear resistance and tool life, owing to a coarse grain size 
and smaller cobalt content in the surface region than in the core region with a graded interlayer, 
Ichikawa (2001). Thus the highest hardness was achieved in the surface layer while the toughness 
increased towards the core region. 



Figure 2. Graded ceramic-metal thermal barrier coating (Source: Dollmeier et al. 1995). 


There is no doubt that FGMs offer substantial advantages over traditional composites in the 
design of structural components which require different types of materials to simultaneously satisfy 
various requirements such as load bearing capability, surface protection, and perhaps decorative 
function. In particular, it has been demonstrated that microstructural gradients not only can over- 
come problems associated with traditional composites by smoothing out property discontinuities, 
but can also produce unique and desirable functions, such as focusing light in fibers, channel- 
ing heat in computer chips, and providing biocompatible implant-tissue transitions in biomedical 
engineering applications, for instance. 

However, there are many challenges that must be overcome in order to make full use of the FGM 
potential, Ilschner (1997). Materials researchers agree that the FGM manufacturing capabilities 
and material development are not in a dynamic growth phase comparable to that which occurred 
during the 1960-1980 period for the then-emerging materials such as plastics. The reason is that 
with the present state of knowledge, it is very difficult and expensive to fabricate graded materials 
with various designed functions in large quantities for industrial use. Moreover, at the present 
time, mature design standards and test methods for FGMs are insufficient or unavailable. Just as 
importantly, the multiscale approach required in the analysis, design and optimization of FGMs 
is still not well developed. The present investigation, therefore, was undertaken to address this 
last point, as will be discussed in more detail later. The literature review of the following section 
provides the motivation for the undertaken study. 
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1.1 Literature Review 


In this section, a brief overview of the FGM development is given for completeness, followed by 
a discussion of theoretical modeling methods used to analyze FGMs. The emphasis is put on the 
latter in view of the focus of this investigation which addresses one aspect of multiscale modeling 
of FGMs, as discussed in the next section. 

With the realization of the important role that FGMs will play in many fields, many researchers 
have devoted a substantial amount of time and energy in recent years to the broad area of FGMs, 
ranging from the different fabrication techniques and characterization methods to the development 
of theoretical models for the response of these materials. The development and exploration of 
FGMs was initiated in Japan under a national project during the period 1987-1991. The aim of 
this government-sponsored project was to develop FGMs for future aerospace applications. Fol- 
lowing this, similar projects and forums for FGM research and development were initiated in other 
countries, in particular in Germany, Finland, China, the United States and Russia. Research on 
functionally graded materials around the world increased rapidly after the First International Sym- 
posium devoted solely to these materials, which was held in Sendai, Japan in 1990. The second 
International Symposium on FGMs was held in San Francisco in 1992, followed by the Third Inter- 
national Symposium held in Lausanne, Switzerland in 1994. These symposia, which were organized 
by the International Forum on FGMs, are now held every two years on a regular basis with the 
next one scheduled for the fall of 2002 in Beijing, China. 

As a result of the above activities, the overall research and development of FGMs has expanded 
to many areas. Consequently, many papers on different aspects of FGMs have been published in 
standard journals, special issues of journals devoted solely to FGMs, conference proceedings, and 
monographs, cf. Pindera et al. (1994a, 1995a, 1997), Needleman and Suresh (1996), Yamanouchi et 
al. (1990), Ilschner and Cherradi (1994), Shiota and Miyamoto (1997), and Suresh and Mortensen 
(1998). Therefore, a comprehensive review of the different research activities is outside the scope 
of this report. Therefore, in keeping with the investigation’s focus, an overview is provided below 
of the different approaches employed by researchers to model the thermomechanical response of 
FGMs. 

The review of the various techniques employed to analyze the response of FGMs under thermal 
and mechanical loading reveals two distinct and fundamentally different computational strategies. 
These have been termed uncoupled and coupled micro-macromechanical approaches by Aboudi 
et al. (1993). The uncoupled approaches do not explicitly couple the material’s heterogeneous 
microstructure with the structural global analysis. Local effective or macroscopic properties at a 
given point within the FGM are first obtained through homogenization based on a chosen micro- 
mechanics scheme, and subsequently used in a global thermomechanical analysis, Figure 3. The 
response of an FGM under specified loading at a homogenized material point with equivalent prop- 
erties is obtained by solving the governing differential equations (partial or ordinary) with variable 
coefficients which represent the spatially macroscopic material property variation obtained from 
homogenization. Sometimes, the spatial property variation is represented as piece-wise uniform 
which is accomplished by discretizing the microstructure into layers or strips with uniform proper- 
ties that change gradually from layer to layer. Clearly, by employing only effective or homogenized 
properties in the global FGM analysis, the microstructural details can only affect the FGM in an 
effective or average way. It is in this sense that the approach is uncoupled. 

The key assumption made in the uncoupled approach is that the heterogeneous microstructure 
of an FGM can be replaced by an equivalent continuum with a set of macroscopic properties that 
vary with spatial coordinates in a manner commensurate with the material’s heterogeneity. The 
various micromechanical approaches that may be used to calculate homogenized material properties 
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Figure 3. Homogenization-based analysis of FGMs. 


include Reuss and Voigt hypotheses, self-consistent schemes, Mori-Tanaka method, concentric cylin- 
der models, and the Generalized Method of Cells, as discussed by Aboudi (1991). Therefore, the 
use of the uncoupled micro-macromechanical approach in analyzing the response of FGMs is im- 
plicitly based on the assumption of the applicability of the representative volume element (RVE) 
concept. However, the definition of an RVE remains to be established theoretically for materials 
with spatially variable heterogeneous microstructures. In addition, ignoring the possibility of cou- 
pling between local and global effects itself may result in potentially erroneous predictions. This 
is especially relevant when the thermal and mechanical field variable gradients are relatively large 
compared with the characteristic dimension of the inclusion phase or the number of inclusions is 
small, or the dimension of the inclusion phase is large relative to the global dimension of the com- 
posite, Aboudi et al. (1993, 1994). Therefore, this type of approach should be used in limited 
circumstances such as when the size of the inclusion is relatively small compared with the overall 
size of composite or the number of the inclusions is large, Pindera et al. (1995b). 

The alternative to the uncoupled approach is the coupled approach wherein the micro-macrostructural 
interaction is explicitly taken into account. This implies that the local micromechanics analysis 
must be combined with the global or structural analysis in a way that connects the two scales 
directly. One way to do this is to solve the entire problem directly by accounting for the presence of 
every local heterogeneity. In the presence of many inclusions, the solution of the given boundary- 
value problem in an exact sense is prohibitive and often impossible, necessitating development of 
approximate computational strategies. 
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1.1.1 The uncoupled macrostructural approach 

The uncoupled approach is the simpler of the two approaches in the analytical sense since it 
treats the micro and macrostructural problems separately. As a result, many particular types 
of problems involving FGM applications have been solved using this approach. One of the most 
important applications involves thermal barrier coating (TBC) problems. As discussed in the 
preceding section, internal stresses can be reduced and fracture toughness can be enhanced by 
appropriate choice of microstructural gradation. For instance, Fukumoto et al. (1994) addressed 
the influence of coating configuration on cyclic thermal shock fracture behavior of plasma sprayed 
coatings. An extensive review of the different approaches employed to evaluate temperature and 
stress fields in TBCs has been recently provided by Noda (1999). In the presence of inelastic 
effects, the approaches range from one-dimensional finite-difference analyses, Kokini and Choules 
(1995), to two-dimensional finite-element analyses based on layer- wise discretization of the coating 
microstructure with piece- wise uniform properties calculated using simple rule-of-mixtures or Mori- 
Tanaka estimates, such as those by Jian et al. (1995) or Delfosse et al. (1997). Elastic analyses 
often use continuously varying properties of particular functional form, thereby facilitating either 
closed-form analytical solutions, or numerical solutions based on the Green’s function approach, 
such as those presented by Nomura and Sheahen (1997) and Sutradhar et al. (2002). 

Another important application of FGMs involves joining structural components made of dis- 
similar materials using graded transition regions to reduce the large interlaminar stresses that arise 
at the free edge due to the material property mismatch. Analyses of these problems have been 
conducted by Drake et al. (1993) and Williamson et al. (1993) using the finite-element approach 
to study the effect of material gradation on the free-edge peel stresses in layered cylinders with 
gradual material variation achieved by using slightly different material properties in each cylindri- 
cal layer. Plasticity was demonstrated to be an important factor in influencing the free-edge stress 
fields. Suresh et al. (1994) demonstrated similar effects in analyzing the response of elastoplastic 
biomaterial strips subjected to cyclic thermal loading. 

Fracture analysis is also an important area of FGMs. In a sequence of papers, Erdogen and co- 
workers extended the elasticity techniques for crack- analysis problems of inhomogeneous materials 
to functionally graded materials, Delale and Erdogan (1983), Erdogan (1985), Erdogan et al. (1991). 
The focus of these studies were the different failure modes that occur in graded TBCs which are 
initiated by free-edge cracks, and cracks parallel or perpendicular to the coating surface, as discussed 
by Kokini et al. (1996) or Kawasaki and Watanabe (1997). These were subsequently analyzed by 
Erdogan (1995), Erdogan and Wu (1996), Lee and Erdogan (1998), and most recently by Jin and 
Paulino (2001). 

Other applications of FGMs include the important problem of identifying optimum microstruc- 
tural designs which produce the best temperature and stress distributions under given loading 
conditions. Various analytical solutions have been developed for specific boundary-value problems 
(typically involving axisymmetric fields and out-of-plane shearing) which can be incorporated into 
optimization algorithms for design purposes. These analytical solutions are also useful as bench- 
mark solutions for the validation of finite-element, finite-difference or boundary-element techniques. 
For instance, Horgan and Chan (1999) developed elastic solutions for problems involving pressur- 
ized hollow cylinders, rotating disks, and bars under torsion with continuously graded isotropic 
constituents. In order to enable fiber-matrix interface/interphase design and optimization in metal 
matrix composites, Pindera et al. (1994b) extended the previously-developed thermoplastic solu- 
tion to the problem of an arbitrary layered concentric cylinder under combined thermomechanical 
axisymmetric loading, Pindera et al. (1993), by incorporating the capability to treat the interfacial 
layers as two-phase composite materials. The results provided insight into designing/engineering 


NAS A/CR— 2002-2 11910 


6 



homogenized and functionally graded interfaces in advanced metal-matrix composites with the ob- 
jective of generating preferred fabrication-induced residual stress distributions. Salzar and Barton 
(1994) subsequently incorporated this solution strategy into a commercial optimization code and 
demonstrated the efficiency with which optimum interfacial architectures could be identified. 

1.1.2 The coupled micro- macrostructural approach 

The coupled micro-macromechanical approach circumvents the limitations of the standard uncou- 
pled approach in FGM applications by explicitly coupling the local (microstructural) and global 
(macrostructural) effects. The coupled analysis can be carried out using the finite-element tech- 
nique. However, this can be prohibitively expensive due to the refined volume discretization required 
for convergence when dealing with even simple FGM microstructures. For this reason, alternative 
approximate computational techniques are needed which are sufficiently accurate to compute the 
micro-macrostructural interaction without requiring excessively large storge. 

The representative of such a coupled micro-macromechanical approach is the higher-order theory 
for functionally graded materials (HOTFGM) developed by Aboudi and co-workers in a sequence of 
papers and reports originating in 1993, Aboudi et al. (1993). Summaries of the different stages of 
the higher-order theory have been provided by Pindera et al. (1995c, 1998), and most recently by 
Aboudi et al. (1999). The theoretical framework is based on volumetric averaging of the various field 
quantities, satisfaction of the field equations in a volumetric sense, and imposition in an average 
sense of boundary and interfacial conditions between the subvolumes used to characterize the 
composite’s functionally graded microstructure. The theory has been validated through comparison 
with the results obtained from finite-element analysis, Pindera and Dunn (1997), and subsequently 
applied to the following technologically important problems: 

• Investigation of the effect of microstructure on thermal and stress fields in MMC plates and 
cylinders 

• Investigation of the use of functionally graded architectures in reducing edge effects in MMC 
plates 

• Optimization of functionally graded microstructures in MMC plates and cylinders 

• Development of guidelines for the design of special coatings in exhaust nozzle applications 

• Investigation of the microstructural effects in functionally graded TBCs 

The results obtained so far have demonstrated that this theory is an accurate and easily imple- 
mentable tool in the analysis and design of FGMs. Furthermore, comparison of the results obtained 
from the standard micromechanics approach with those of HOTFGM has demonstrated the need 
for a theory like HOTFGM, which explicitly couples the micro (local) and macro (global) effects 
within a unified analysis. 

The recent developments of HOTFGM include extension to cylindrical coordinates to enable 
analysis, design, and optimization of structural components found in aircraft engine applications, 
Pindera and Aboudi (2000). 

1.2 Objectives of the Investigation and Outline of the Report 

Despite the demonstrated utility of the higher-order theory in applications where the micro- 
macrostructural coupling effects are important, the theory’s full potential has not yet been attained. 
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This is because the theory as originally formulated is computationally intensive. This, in turn, lim- 
its the types of problems that can be solved due to the large number of equations required to mimic 
realistic material microstructures. Therefore, the objective of this investigation was to construct 
a basis for an efficient reformulation of the higher-order theory, and subsequently implement it in 
order to reduce the number of equations required in the solution of a given thermal boundary- value 
problem. The basis for the efficient reformulation is the demonstration that the same results are 
obtained when the volume discretization used in the original version is simplified by eliminating 
the need to differentiate between generic cells and subcells. This is called a user-friendly formu- 
lation which, in turn, provides the basis for an efficient reformulation based on the local/global 
conductivity matrix approach. The efficient reformulation admits spatially variable thermal con- 
ductivity within individual subcells, thereby making it possible to treat microstructures with two 
fundamentally different length scales. 

The theoretical framework of the original two-dimensional formulation of HOTFGM, limited to 
thermal boundary-value problems, is briefly described in the first part of Section 2 , with the sec- 
ond part devoted to the development of the user-friendly formulation. Based on the demonstrated 
feasibility of the simplified user-friendly formulation, an efficient reformulation of HOTFGM with 
variable thermal conductivity is developed in Section 3 using the local/global conductivity matrix 
approach, which is an extension of the local/global stiffness matrix approach described by Pindera 
(1991) in the context of purely mechanical boundary-value problems. Analytical solutions used 
to validate both the user-friendly and efficient reformulations are developed in Section 4. Volume 
discretization sensitivity and validation studies, as well as a practical application of the developed 
efficient reformulation, are provided in Section 5. Section 6 summarizes the present accomplish- 
ments and discusses future desirable extensions of the developed computational capability. 


2 Higher- Order Theory: Thermal Problem 

2.1 Outline of the Original Formulation 

The original formulation of the two-directional higher-order theory for functionally graded ma- 
terials (HOTFGM-2D) was outlined by Aboudi et al. (1996). The basic idea was to develop a 
two-dimensional framework for the response of composites which are functionally graded in two 
directions, accounting for micro-macrostructural coupling effects. As indicated in Figure 4, the 
heterogeneous composite model has finite length L and width H in the X 2 — x% plane and extends 
to infinity in the x\ direction. The thermal loading applied on the boundary can be specified in 
terms of temperature or heat flux distribution. The microstructure of the heterogeneous composite 
is discretized into N q and N r generic cells (g, r) in the intervals 0 ^ x^ ^ H, 0 ^ x% ^ L in the 
X 2 ~ plane. The indices q and r represent the location of the generic cell in the X 2 — X 3 plane and 
they can vary from q = 1,2... to N q and r = 1,2,... to N r . Each generic cell is further discretized 
into four subcells denoted by the indices (/?y), Figure 4, where the range of each index /?, 7 is from 
1 to 2. These indices indicate the relative position of a given subcell within the generic cell. The 
dimension of the generic cell along the x\ direction is infinite whereas the dimensions along the X 2 

(q) (q) (r) (r) 

and 3 axes, Jq, hq and l\ J ' 2 , respectively, can vary from cell to cell such that 

N q N r 

H = £ (A («> + ,,<«>) , L = £(,M + ,<->) 

q = 1 r= 1 
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Figure 4. Schematic of the HOTFGM analysis geometry. 


An approximate solution for the temperature field under steady-state thermal loading is con- 
structed on the basis of volumetric averaging of the heat conduction equation, together with the 
imposition of boundary and continuity conditions in an average sense between the subvolumes used 
to characterize the material’s microstructure. The solution to the heat equation is accomplished 
by approximating the temperature field in each subcell of a generic cell using a quadratic expan- 
sion in the local coordinates x^\x^ centered at subcell’s midpoint. The higher-order terms in 
the temperature field expansion account for the local effects caused by the thermal field gradients, 
the microstructure of the composite and the finite dimensions in the functionally graded direc- 
tions, in contrast with the linear expansion used in constructing the micromechanics model for 
periodic multiphase materials called the Generalized Method of Cells (Paley and Aboudi, 1992). 
The higher-order theory is an extension of this micromechanics model to materials without periodic 
microstructures, and employs similar volume discretization and related notation. 

The unknown coefficients associated with each term in the temperature field expansion are ob- 
tained by constructing a system of equations that satisfies the requirements of a standard boundary- 
value problem for the given field variable approximations. That is, the zeroth, first, and second 
moments of the heat conduction equation are satisfied in a volumetric sense. The thermal and 
heat flux interfacial continuity conditions within a given cell, as well as between a given cell and 
adjacent cells, are imposed in an average sense together with the applied boundary conditions. 

The following two subsections outline the thermal problem formulation and its solution based 
on the volume discretization employed in the original higher-order theory in order to prepare the 
stage for a generic cell-free, user-friendly formulation outlined in the second part of this section. A 
detailed exposition was provided by Aboudi et al. (1999). 


2.1.1 Thermal analysis 

Heat conduction equation Assume that the functionally graded material model is subjected 
to steady-state temperature or heat flux distribution on its bounding surfaces in the X 2 — x% plane. 
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Under steady-state heat conduction, the heat flux field in the material occupying subcell (/3j) of 


the ( q , r)th generic cell, in the region 


JP) 


— 2 '73 ’ 


™(7) 


X- 


< ±/ (r) 

— 2^7 ’ 


must satisfy: 


dx ^ dx 


= 0 


(1) 


The components of the heat flux vector in this subcell are obtained from the Fourier’s heat 
conduction law for anisotropic materials, 



_ u ( 07 ) 
K ij 


7) 


(z, j = 2, 3; no sum) 


(2) 


where is the temperature distribution in the subcell (/Jy) of the (g, r)th generic cell measured 

with respect to a reference temperature, and k\^ are the coefficients of thermal conductivity of 

the material in the subcell (/?y), with k\ ^ = k\^ 6 ij (no sum on i) for orthotropic materials. No 
sum is implied by the repeated Greek letters in the above and henceforth. 

The solution to the heat conduction equation in each subcell of every generic cell is obtained 
in conjuction with satisfaction of the continuity and boundary conditions in the manner described 
next. 


Heat flux continuity conditions The continuity of the heat flux vector q(^ 7 ) a t the inter- 
faces separating adjacent subcells within the generic cell (< 7 , r) is fulfilled by imposing the integral 
conditions 
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In addition to the above surface- averaged continuity conditions within the (< 7 , r)th generic cell, the 
heat flux continuity conditions at the interfaces between adjacent cells are given by 


,Z «/2 

/-jW /2 

rhf/ 2 
kh^/2 


,( 27 ) I (Q,r) 


q 2 " l 4 2 ) =4 ,) /2 d *' 3 


Jn) - 


% 


(/32) I (q,r) 




=U 72 


dx<7 = 


U7/2 
/-4V 2 
rhf/ 2 

l-h^/2 


% 


(I7) i(5+l,r) 


l 4 1) =-/i^ +1) /2 


dx 


~(7) 




0 31 ) |(«X+ 1 ) 


l 4 1 ) =-4 r+1) /2 


dx 


=09) 


( 5 ) 

(6) 


Thermal continuity conditions The thermal continuity conditions at the interfaces separating 
adjacent subcells within the generic cell (< 7 , r) are similar to the corresponding heat flux continuity 
conditions, 


/ 4 r) /2 

/- 4 r) /2 

/•7V 2 


^( 17 ) |(^’ r ) 


a; 2 1)=fe( i 9) /2 


dx7 = 


j-(/3i) I («.»•) dx^ — 

(«),, 2 “ 


'-4 9 72 


r4 r) /2 

/- 4 r) /2 

/•7V 2 

l-h^/2 


^(27) |(5> r ) 


( 2 )__,(<z) 


=-h^/2 


dx 


*( 7 ) 


j7(/32) |(<?70 


«( 2 )_ 


=- 4 r) /2 


dx 


09) 


( 7 ) 

(8) 
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Similarly, the thermal continuity conditions at the interfaces between neighboring cells are given 
by 


4 r) / 2 


7(7 2 

r44 / 2 


^( 27 ) |(?5' r ) 




=W 72 


= 


4 r) /2 


-4^/2 


j^(/32) ifer) 


/-7 } /2 


J>(l7) |(« , + 1 > r ) 


s ( 1) =-/ l (‘ J+1) /2 


da: 


(7) 


^ ( 2 ) = 4 r) /2 


da;)/ 3 -* = 


* / 2 t(^) | ( ^ +1) (r+1)/ dxf 
(<z ) /0 '4 )= -4 + / 2 2 


7 -^72 


(9) 

(10) 


Boundary conditions The solution for the temperature field must satisfy the boundary condi- 
tions at top and bottom, left and right surfaces of the graded material, i.e., on the external surfaces 
of the boundary cells q = 1,A^; and r = l,iV r . When q = 1 and q = N q , the temperature in 
the cells (l,r) and (N q ,r) at the X 2 = 0 and X 2 = H (i.e., left and right) surfaces, respectively, 
must equal the applied temperature. As in the case of the interfacial continuity conditions, these 
boundary conditions are imposed in an average sense, 


^7/2 

-472 


y(i7) |( 1 > r ) 


*(i)_ 


— _ i/T 
— 2 n l 


( 1 ) 


dx 


*( 7 ) _ 


r4 r) /2 

/- 4 r) /2 


' <7 * d;rl 7 * 


ie/t 


(ii) 


r4 r) /2 

/-4 r) /2 


^( 27 ) |(-^b, r ) 


( 2 )_ , I..W 7 ) 
— ^ 2 


dx 


(7) _ 


r4 r) /2 

-472 


"'■ 7 l da’ 7 


right 


( 12 ) 


where r = 1, , 7V r . 

Similarly, when r = 1 and r = 7V r? the temperature in the cell (< 7 , 1 ) and (< 7 , A^ r ) at the £3 = 0 
and X 3 = L (top and bottom) surfaces, respectively, must equal the applied temperature, 


-44/2 


/- 7 9) /2 


j-(/ 3 l) 1(3,1) dr^ — 

7( 1 )__ i/C 1 ) “ X 2 — 

X 3 — 2 1 1 


rhf/ 2 

l-h ^/2 




(13) 



rn(/3 2 ) 

1 1 ^( 2 ) , l/(Mr) 

x 3 — ^ 2 6 2 




rp{l) 
- 1 top 



(14) 


where q = 1 , , TV^. 

Alternatively, we can also apply mixed-boundary conditions involving temperature and heat 
flux on different portions of the boundary. 


2.1.2 Solution of the thermal problem 

The temperature distribution in the subcell (/?y) of the (< 7 , r)th generic cell, measured with respect 
to a reference temperature T re y, is denoted by T^ 7 ). The temperature field is approximated by a 
second-order expansion in the local coordinates x^\x^as follows (omitting the indices (< 7 , r) for 
notational simplicity): 


'T'iPl) — T 7 ^ 7 ) 1 ~(P) r r(P'y) 1 ~{l)rn{Pl) 

1 1 ( 00 ) 17 x 2 ^(10) 77 o ^(01) 


+ i(3x 


2 


7 ) 2 


\rp{Pl) 
~ ' ( 20 ) 


+ 7 


^(7 ) 2 

3 


,b ) 2 


' ( 02 ) 


(15) 


In the two-directional version of HOTFGM, the temperature field varies only in the X 2 — X 3 plane, so 
the temperature field expansion terms in the x\ direction are not included. In the above temperature 
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field approximation, we can see clearly that there are five unknowns, where is the volume- 

averaged temperature within the subcell, and t|^^ (m, n = 0 , 1 , 2 with m + n < 2 ) are higher-order 
terms which are determined from the conditions outlined previously. The reason for the absence of 
the product terms (i ^ j) is due to the averaging procedure employed in the higher-order 

theory. 

In summary, in each generic cell there are four subcells and in every subcell the five unknown 
coefficients, 

[T m > 7(io) ^ 7(oi) , 7(20) 5 T ( o2)] (/37) 

describe the temperature field. Therefore, there are altogether 20N q N r unknown quantities which 
must be determined for a composite with N q rows and N r columns of cells containing arbitrarily 
specified materials within their subcells. These unknown coefficients are determined by first sat- 
isfying the heat conduction equation, as well as the first and second moment of this equation in 
each subcell in a volumetric sense. Then the temperature and heat flux continuity conditions are 
applied in an average sense at the interfaces separating adjacent subcells as well as adjacent cells. 
Finally, the boundary conditions are imposed in an average sense. Following this procedure, we 
obtain 20N q N r equations for the 20N q N r unknown coefficients in the temperature field expansion 
as summarized below. 


2.1.3 Heat conduction equation 

In the course of satisfying the steady-state heat equation in a volumetric sense, the heat flux 
quantities given below are defined for convenience 


n ^ = 


1 

4 (?r) 

A m 


+ 9 + r+& 


03) 


P) , 


( 7 ) ' 


l-h (<?) /2 /2 

n m /z \i) /z 


[x^ ) m (x^ ) n q[^ dx dx^ 


(16) 


where m, n = 0 , 1 , or 2 with m + n < 2 , and is the cross-sectional area of the subcell 

(/?7) in the (< 7 , r)th generic cell. For example, when m = n = 0, Q^qo) i s th e volume- averaged 

value of the heat flux component qf 7 ^ in the subcell, while for other values of (m, n), higher-order 
volume integrals of the heat flux components are obtained. These flux quantities can be evaluated 
explicitly in terms of the coefficients T^ 7 ^ by performing the required volume integration using 
Eqs. (2) and (15) in Eq. (16). This yields the following non- vanishing zeroth-order and first-order 
heat fluxes expressed in terms of the unknown coefficients in the temperature field expansion. 


Q 

Q 

Q 

Q 


(07) _ 

2(00) 

u (07)71(07) 
"'2 1 (10) 

(17) 

(07) _ 

2(10) 

A q ) 2 

7,(0 7) 0 rp(P 7) 

2 4 (20) 

(18) 

(07) 

3(00) 

1.(07) 7^(07) 

o 1 (01) 

(19) 

(07) _ 

3(01) 

»(r)2 

7,(07) 7 t.(07) 

K 3 4 ^ (02) 

( 20 ) 


Satisfaction of the zeroth, first and second moment of the steady-state heat equation, Eq. ( 1 ), in 


volumetric sense results in the following four relationships among the first-order heat fluxes Q 


(£ 7 ) 
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in the different subcells (/?y) of the (g, r)th generic cell, where (/?y) assumes all permutations of 
the integers 1 and 2, 


[Qzm/hp + Qfm 


/ifi q ' r) = 0 


( 21 ) 


Heat flux continuity equations The equations that ensure heat flux continuity at the subcell 
interfaces, as well as between individual cells, associated with the X 2 and directions, Eqs. (3)- (6) 
are given by 

[-12Q^» 0 )/A, + Og> 0) - 6Q%/fc2] M - [«2(lo) + <(lVM h “ 1,r| = 0 (22) 

^ 2 ( 0 ( 0 ) + 2 ^ 2 ( 0 , 0 ) “ 3 < ? 2 ( T , 0 ) + 2 ^ 2 ( 0 , 0 ) + 6 < ? 2 ( l , 0 )/ / i ' 2 l (9 ~ 1 ’’' 1 ^ 0 ( 23 ) 

+ «S$» - * q%\M w) - C) + = 0 < 24 > 

l-Qfm + = 0 < 25 > 


Equations (22)- (25) provide eight relations among the zeroth-order and first-order heat fluxes. It 
is observed that Eqs. (22)-(25) together with Eq. (21) can be expressed in terms of T mn^ by using 
Eqs. (17)-(20). 

Thermal continuity equations The equations that ensure the continuity of temperature at the 
interfaces between the subcells, as well as between the neighboring cells in the X 2 and xs directions, 
obtained from Eqs. (7)-(8) and Eqs. (9)-(10), are as follows: 


l T m + s'*! 2 ’™’ + -A = pffi? - (VS + i h 2 T Koll i, ’ r> 


( 10 ) 


’( 20 ) 


( 00 ) 


( 10 ) 


( 20 ). 


PfS! + = pr$> - 

\rri(01) . I 7 rri([31) , ^ i2rri([31) 1 [q r) |>ri(/32) ^ j rp((32) 1 i2rp((32)^ (q r ) 

Hrocn + 0^(01) + 7 / i 7 rn 2 ^ r ~ Ton) “ o^ron + ih 1 ^ \ 


-( 00 ) 


( 02 ) 


( 00 ) 


L (01) 


( 02 ) 


[ rp\P^) I ^yn(/32) ^ j2 rp{(32) i (g,r) rrp({31) __ ^7 rp{{31) , ^ ?2 ^(/^l) i (g ?r +l) 

^ 2^ 2 ( 01 ) ^ 4^ 2 (02) J — (oo) 2^ (01) ^ 4^ (02) 1 


Equations (26)- (29) provide us with the remaining eight relations governing the temperature 
within interior generic cells. 


(26) 

(27) 

(28) 
(29) 
field 


Governing equations The heat conduction equation, Eq. (21), together with the heat flux 
and thermal continuity conditions, Eqs. (22)- (25) and (26)- (29), respectively, form 20N q N r linear 
algebraic equations which govern the 20^A^ r microvariables Tmn^ in the four subcells (/?y) of an 
interior cell (g, r); q = 2, ...N q — 1, r = 2, ...N r — 1. For the boundary cells q = 1, Eq. (21) as well 
as Eqs. (24)- (29) are still applicable. However, the heat flux continuity conditions between a given 
cell and the preceding cell given by Eqs. (22) and (23) are not applicable. They are replaced by 
the applied temperature or heat flux conditions at the left boundary X 2 = 0, i.e., Eq. (11), and 
the condition that the heat flux at the interface between subcell (I7) and (2y) of the cell (l,r) is 
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continuous. For the boundary cells q = N q , the previous equations are applicable except for Eq. 
(27), which is replaced by the specific temperature applied at the right boundary X 2 = 77, be., Eq. 
(12). Similar situation holds for the boundary cells r = 1 and r = N r . 

The governing equations for the interior and boundary cells form a system of 20N q N r algebraic 
equations in the unknown coefficients T The solution of these equations determines the temper- 
ature distribution within the functionally graded composite subjected to the boundary conditions 
specified by Eqs. (11)-(14). The final form of this system of equations is symbolically represented 
below: 

kT=t (30) 


where the structural thermal conductivity matrix k contains information on the geometry and 
thermal conductivities of the individual subcells (/3y) in the N q N r generic cells (g,r). The thermal 
coefficient vector T contains the unknown coefficients that describe the temperature fields in each 
subcell. It is given by 


T = [T 


(ii) 

ii, 


T 


( 22 ) ] 
N q N r J 


where 



Finally, the thermal force vector t contains information on the boundary conditions. 


2.2 User-Friendly Formulation of HOTFGM-2D 

The motivation for the development of user-friendly thermal formulation of HOTFGM-2D is to 
ultimately construct a basis for an efficient reformulation using the local/global conductivity matrix 
approach, which will be outlined in the next section. User-friendly formulation of HOTFGM- 
2D basically follows the solution technique of the original formulation which was outlined in the 
preceding section. The major difference between these two formulations is simplification of the 
equations through the elimination of the generic cell concept which also simplifies the heat ffux 
continuity conditions. 

In the user-friendly formulation, therefore, the microstructure of the heterogeneous composite 
is discretized into NpN'j subcells in the intervals 0 < X 2 < 77, 0 < < T, i.e., the X 2 — x% plane, 

Figure 5. The arbitrary subcell (/3y) has length hp, and width l y , where the /? = 1,2, , Np\ 

7=1,2, , Nj. In the original formulation, there are five unknown quantities associated with each 

subcell and four subcells within each generic cell, therefore, 20N q N r unknown coefficients must be 
determined by establishing a 20^A^ r system of equations involving field equations and continuity 
conditions, together with the imposed thermal or heat flux boundary conditions. In the user-friendly 
formulation, there are still five unknown quantities in the temperature field expansion associated 
with every subcell, because the same temperature expansion is employed. Since the total number 
of subcells is NpN 1 ^ 5 NpN 1 equations must be constructed for the determination of the unknown 
quantities. Although the number of unknowns in each subcell and the total number of equations 
remain the same in the user-friendly formulation (since 2 N q = Np and 2N r = TV 7 ), the simplified 
notation leads to a more systematic method of assembling the final system of equations for the 
determination of the unknown coefficients in the temperature field expansion. More importantly, 
it also provides a basis for the efficient reformulation described in the next section. 
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2.2.1 Thermal analysis 

As before, the temperature distribution in the subcell (/3y) is approximated by 


r r((3'y) 7 ) 1 ^{(3)rri((3 7 ) . ^{l)rp{0l) , }_( 0^,(0) % 

± — - L fnn\ "T ^9 -Wirrt "T" -t "T" 2 V ox 2 


-( 00 ) 


( 10 ) 


-( 01 ) 


h 2 

0_\rn(0l) 

4 ; (20) 


+ o( 3 4 7)2 


Z\rp(Pl) 

4 j ( 02 ) 


(31) 


Heat conduction equation As in the original formulation of HOTFGM- 2 D, under steady-state 
he heat flux field in the material occupying subcell (^ 7 ) in the region \x\\ < 00 , 
< ^/ 7 must satisfy Eq. (1) in a volumetric sense, 


heat conduction, tl 

=09) 

x 2 

— 2^M 

=( 7 ) 

x 3 


M/ 37 ) <9^, 


0 M 7 ) 0 M 7 ) 

'^2 + %^ )dV m =0 


dp) 


dx 


=(7) 


(32) 


The heat flux components and q^ 1 ^ can be expressed in terms of the unknown coefficients 
in the temperature field expansion, Eg. (31), using Eq. (2), 


«r> = **r> 


QT^Pi) 


dx . 


( 0 ) 


dT^ 


dx . 


M) 


_b,{0l) (rn{0l) I rnifi'y) \ 

— ^2 ( 1 ( 10 ) ° x 2 1 ( 20 ) / 

_ _iAPt) (rp(Pl) 1 o^l)rn{0l)\ 

— ^3 (1(01) “T" ox 3 1(02)/ 


(33) 

(34) 
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Therefore, substituting Eqs. (33) and (34) into the steady-state heat conduction equation, Eq. 
(32), yields, 

rn(fil) , u(Pl)rji{Pl) Q /OC\ 

^2 1 ( 20 ) ^ ^3 1 ( 02 ) — U v 00 ) 

This is the same equation as Eq.(21) in the original formulation at the subcell level, after replacing 
the temperature microvariables with the volume- averaged heat fluxes using Eqs. (18) and (20). 


Heat flux continuity conditions between adjacent subcells Using the new generic cell-free 
notation, the heat flux continuity condition between adjacent subcell’s along the X 2 direction is 



(A7) 


X 2 — hp 1 2 




(0+1,7) 

<?2 


\x 2 =—hf 3 +1 /2 



(36) 


Substituting for the heat flux components in the x 2 direction using Eq. (33), the above equation 
becomes: 


jJ/3,7) (rp(P,l)i _|_ 3hp T (0,-y) , , _ (/ 3 +l, 7 ) / t (/3+1,7), _ 3hf3+ lm^+ 1 , 7 ). \ 

^2 I” 1 (10) l 7 7“ 2 "A 2 *)) ‘7/ — "'2 1 J (10) l 7 2 "A 2 *)) ( 7/ 


Similarly, the heat flux continuity condition along the X 3 direction is 


rhp/2 


% 


(Pn) 


\X 3 =l ^/2 


dxi^ = 


J-hp/2 

which, upon using Eq. (34), becomes 


r h p / 2 

Z-^/2 


(0,7+!) | ^( 0 ) 

^3 l® 3 =-Z 7+ i/2 ax 2 


(''W + ! % - %) 


(37) 


(38) 


(39) 


Thermal continuity conditions between adjacent subcells Using the new notation, the 
thermal continuity condition at adjacent subcells interfaces along the X 2 direction is, 



/A / 2 
'~h/ 2 


U.,-^,/2 dx { p 


(40) 


Upon substituting for the temperature field defined by Eq. (31) and performing the required 
integration, the above condition reduces to 


77 ( 0 7)7 | ^ 0 ^ 77 ( 07 ) | ^0^7 t (0 7 ) _ ( 0 + 1 , 7 ), ^ 0 + 1 ^ 7 77 ( 0 + 1 , 7 ) ^0+1^7 ^,( 0 + 1 , 7 ) //j-ix 

1 ( 00 ) 1 ' 2 ( 10 ) ^ 4 ( 20 ) ( 00 ) / 2 ( 10 ) ^ 4 ( 20 ) 

Similarly, in the £3 direction, the thermal continuity condition is, 

[ h0/2 \x 3=h/2 dxf = [ hp/ 2 T^’7 +1 ) |* 3= _ W2 ctef (42) 

4 — hp/2 4 — hf3/2 

which, upon using Eq. (31), becomes 


rp(Pl ) h i W}P_rp{P^) , W^0_rp{0 7 ) _ T (/3,7+l) L _ ^ +1 7? rpjP^+1) , 7+j_^ ^(/VH-l) 

J ( 00 ) ' t / 3 7" 2 J ( 01 ) 7" ^ J ( 02 ) — J ( 00 ) 'T 3 2 J ( 01 ) 7" 4 j ( 02 ) 


(43) 
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Boundary conditions The temperature or heat flux at the exterior surfaces of the boundary 
subcells (/?, 1) and (/?, TV 7 ) at X 3 = 0 and X 3 = L, Figure 6 , respectively, must equal the applied 
boundary condition (given temperature or given heat flux) in an average sense, 



rhf 3/2 

/ T m) 

J-hp/2 


\X3 = — l±/2 



rn((3, 1) 
i (00) 


_ (p04); 

1 (01) 11 


/2 + T<* 1 >(?/4 = 



n~i{P) 

1 bot 



1 bot 


(44) 


or 




X3=—li/2 


dxi® = — h 


04 ) frp(P,l) 

O(oi) 






(45) 


and 


f hf>/2 T 0 ,iv 7 ) 

J-hp/2 


X3=In~,/2 



rp{P>N-y) 

1 m 


rp(/3,N^)j 
+ i (01) <JV 7 


/2 + T, 


(/W,2 

( 02 ) 1 


at 7 / 4 


or 



rpifi) 
- 1 top 



rn{P) 
- 1 top 


(46) 



_(/W 

% 


® 3 =£ jv 7 /2 



,(/3,lV 7 ), (/3,1V 7 ) 

K 3 O(oi) 


+ > 


r^/3/2 

/ 2 )= / 

J-h p /2 


iw dx 


i 0 d ( 47 ) 


Similarly the temperature or heat flux in the boundary subcells ( 1 , 7 ) and (TV#, 7 ) at X 2 = 0 
and X 2 = H must be equal to the applied boundary condition (given temperature or given heat 
flux) in an average sense, 



X 2=— hi /2 



'7~'(b7) 

( 00 ) 



V2 + T ( ( 2 ^/4 



T (7) 

1 left 



7) 

1 left 


(48) 


or 
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(!,7) 

Q 2 


X 2 =—hi /2 



( 1 >7)/rp(l,7) 
( 10 ) 




:MmT, 


(1,7) 


( 20 ) 



Qteft* 1 * s 7) 



(49) 


and 


rl'j/Z 

I-IJ2 


T^'% 2=hNp/ 2dx™ 


t (N, 3,7) 
J (00) 


+ T^ ) h N J2+T' : ^hlJi 


L (20) 



rp{ l) 

± right 



= T 


(7) 


right 

(50) 


or 



„W3,7) 

^2 




/2 


^(7) 


dx, 


tA n Pi 7) (rp(.Npi7) 
k 2 K 1 ( 10 ) 


+ 3/iAT>27 


(JV/3,7) 

( 20 ) 





bright 


bright ( 51 ) 


2.2.2 Assembly of the governing equations 

As stated before, there are five unknown quantities in the temperature field expansion associated 
with an arbitrary subcell (/3j). They can be represented by the vector defined as follows, 

T (/37) = 7(00), T m , T (01) , T (20) , T m }^ (52) 

Therefore, we need to construct five equations for every subcell to solve for the unknown coefficients. 
Two equations involve temperature continuity in the X 2 and X 3 directions, another two involve heat 
flux continuity in the X 2 and x% directions, and the remaining equation is the steady-state heat 
conduction equation. These will be arranged into an organized system of equations along each row 
sequence and column sequence of subcells representing the functionally graded material. 

Let the assembly of the system of equations proceed first from left to right (in the X 2 direction) 
along each row of subcells starting from the first row, Figure 7. Subsequently, we proceed from 
bottom to top (in the x% direction) along each column of subcells starting from the first column. 
Finally, we take care of the steady-state heat equation within each subcell along either a row 
sequence or a column sequence. In summary, we obtain two equations at each interior interface 
along either direction after rearrangement, one involving temperature continuity and the other 
associated with heat flux continuity. We also obtain one steady-state heat equation for each subcell 
along either a row sequence or a column sequence for a total of five equations. 



Figure 7. Governing equation assembly. 
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Continuity and boundary conditions along x<i direction For a given 7, the left boundary 
condition, Eq. (48) or (49), can be rewritten in matrix form as follows, 


K b 


rp(i 7 ) 


— or 

“ 1 left ° r 9/ e /t 


(53) 


where depends on the applied boundary condition, and T^j t and q^J t are the surface integrals 
of the temperature and heat flux on the external surface of the boundary subcell ( 1 , 7 ) given by 
the right side of Eqs. (48) and (49). If the applied boundary condition is temperature, then 

Kf 7 = [1 - hi/2 0 h\/A 0] 

Alternatively, if the applied boundary condition is heat flux, then 

Kf 7 = [0 - 4 1,7) ° - 3/2 4 1,7 A 0 

The definition of T^i) [ s given by Eq. (52). 

Proceeding into the interior, we can rewrite the temperature and heat flux continuity conditions 
at the interface separating two adjacent subcells (/?, 7) and ( j 3 + 1, 7) in the interior, Eqs. (37) and 
(41), along the X 2 direction as follows, 


K r+ 

7 

\C T ~ 

)0+l,7 




' 0 ' 

k h + 

tsH 

^£+1,7 


'X’(7 3+1 >7) 


0 


(54) 


where 


K r+ - 
^£7 “ 


K 


T- 

'£+1,7 


1 hp/2 0 h 2 ^/ A 0 

-1 hp+ 1/2 0 — h 2 p +1 /A 0 


K H + — 
^£7 “ 


0 -k^ l] 0 — 3/2hpk 2 


(Pi) 


and (5 = 1 , 2 ,, 


K 


H- 

73 + 1,7 


0 A:f +1 ’ 7) 0 - 2,/2hp +1 kf +la) 0 


,N« 


1 . 


Finally, the right boundary condition, Eq. (50) and (51), can be rewritten in matrix form as, 


7 


f p(^ V /3’7) 


_ rp{ 7) 

± right 


-(7) 

OT tight 


(55) 


where K depends on the applied boundary condition, and T riy f d and (i riy f d 


iljht and bright are the surface 
integrals of the temperature and heat flux on the external surface of the boundary subcell (TV#, 7 ) 
given by the right side of Eqs. (50) and (51). If the boundary condition is temperature, then 


K]V^7 = 


1 h N J 2 0 h 2 N J 4 0 


Alternatively, if the applied boundary condition is heat flux, then 


K 


N 0 ,7 


4^’ 7) 0 - 3 / 2 k { 2 N ^ j) h N 0 
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Therefore, for an arbitrary row 7 , let f3 = 1, Figure 7, and apply the above equations 

from the left boundary to the right boundary, respectively. This procedure produces a system of 

equations defined by the matrix K^, which contains information on the left boundary condition, 
temperature and heat flux continuity conditions at the interfaces between adjacent subcells, and 

the right boundary condition, respectively, along the row sequence. The structure of whose 

size is 2Np x 5 TV#, is shown below, 



K? 7 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

k t+ 

& 

k t_ 

*7+1,7 

0 

0 

0 

0 

0 

k h + 

*P7 

k h ~ 

*7+1.7 

0 

0 

0 

0 

0 

0 

tc t + 

*7+ 1,7 

K r_ 

*7+2,7 

0 

0 

0 

0 

0 

tsH+ 

*7+1,7 

k h ~ 

*7+2,7 

0 

0 

0 

0 

0 

0 

0 



0 

0 

0 

0 

0 

0 

k£ 



Following the same procedure applied to every row, 7 = 1,2, ....TV 7 , we obtain a system of 
equations defined by the matrix A, consisting of the individual matrices that has the form, 

k 7 0000 
0 000 

A = 0 0 Kf 0 0 

0 0 0 0 

0000 

Since the matrix A consists of iV 7 matrices k 7 placed along the main diagonal, its size is 2 NpN^/ x 
5 NpN^. 


Continuity and boundary conditions along X 3 direction Repeating the preceding procedure 
along the X3 direction, the bottom boundary condition, Eq. (44) or (45), for a given /?, can be 
rewritten in matrix form as follows, 




rn{P) 
1 bot 


_(/ 3 ) 

° r %ot 


(56) 


As in the case of the boundary matrix in A, also can be obtained in two ways. When 
temperature is applied to the boundary, then we get 

Kfi =[10 -h/2 0 zJ/4 

Alternatively, when the boundary condition involves heat flux, we get 

Kfj =[ 00 - kf l) 0 2>/2kf : l) h 


and are the surface integrals of the temperature and heat flux on the external surface of 
the boundary subcell (/?, 1) given by the right side of Eqs. (44) and (45). 
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Proceeding into the interior, we can assemble the temperature and heat flux continuity condi- 
tions at the interface separating two adjacent subcells (/?, 7) and (/?, 7+ 1) in the interior, Eqs. (39) 
and (43), along the x% direction as follows, 


where 


K T + 

t>T- 

^, 7+1 


rp(/3 7 ) 


' 0 ' 

K H + 

^37 

i >H- 
^,7+1 


X(/ 3 >7+ 1 ) 


0 


(57) 


77 = [1 0 y 2 0 iy 4] 

Kj- +1 = [-1 0 Z 7+1 / 2 0 -% + J 4] 


- 


0 0 -fcf 7) 0 -3/2 


(^7)7 


K 


H- 

3 , 7+1 


0 0 4 /?,7+1) 0 - 3/2ifef’ 7+1) i 


7+1 


and 7 = 1, 1V 7 — 1. 

Finally, we rewrite the top boundary condition, Eq. (46) or (47), in matrix form as, 



T (/W 


rp{P) 
- 1 top 


or 



(58) 


where T^} and 9^ are the surface integrals of the temperature and heat flux on the external 
surface of the boundary subcell (/?, Ay) given by the right side of Eqs. (46) and (47). If the applied 
boundary condition is temperature defined as T^, then 


K 


B 

/3,iV 7 


1 0 I n J2 0 7/4 


Similarly, if the boundary is subjected to heat flux given by the form q^, then we obtain 


K 


B 

3,1V 7 


0 0 




3/2 kf' N ^l Nl 


Following these steps, for an arbitrary column /?, letting 7 m 1 — » TV 7 , Figure 7, we construct a 
system of equations defined by the matrix which contains information on the bottom boundary, 
interfacial temperature and heat flux continuity conditions, and the top boundary, respectively, 
along the column sequence (i.e., in the X 3 direction). 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

k t + 

•• 737 

t>T- 

" ^3,7+1 

0 

0 

0 

k h+ 

•• -*7*7 

T>ri- 

■■ ^3,7+1 

0 

0 

0 

0 

0 

0 

k t+ 

' 

•• ^3,7+1 

t>T- 

” *&7+2 

T^n- 

•• 73 , 7+2 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

k b 

■■ TsAb ••• J 
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In the above matrix, the designation . indicates (3 — 1 0 matrices, the designation .. indicates 
Np — 1 0 matrices, and the designation ... indicates Np — f3 0 matrices, with the zero matrices 
having the dimensions [1 x 5]. The size of the matrices is 2N~ f x 5NpN~ r 

Extending to all 3,3 = 1. 2, ...., Np, we obtain the matrix B that has the form 


B = 


K 


(7) 


K 


(7) 


K 


(7) 

Np 


with B having the dimensions 2NpN 1 x 51YglV 7 . 


Heat conduction equation in X2 or X3 direction Finally, we write the steady-state heat 
conduction equation, Eq. (35), for an arbitrary (J3 r y) subcell as follows 


where 




= 0 


K*(/?.7) = [0 0 0 - 3 fc? 7) - 3/cf 7) ] 


„(Ay)i 


Assuming that the assembly proceeds from left to right along a row sequence, we obtain the matrix 

*(3) 

K 7 ; with dimensions Np x 5 Np for a given 7 row of subcells 



k* (1,7) 0000 
0 .... 0 0 0 

0 0 k* (/3,7) 0 0 

0 0 0 .... 0 

0 0 0 0 K* (iV/3 ’ 7) 


where (3 = 1 , 2 , Np along the diagonal, respectively. 

Extending the procedure to all 7, we obtain the final system of equations that satisfies the heat 
conduction equation in every subcell. This system of equation is represented by the matrix C, 

which contains the individual matrices K 7 ' along its diagonal as shown below. 


C 


k* (/5) 0 0 0 0 

0 0 0 0 

0 0 k 7 ] 0 0 

0 0 0 0 

0 0 0 0 


where 7 = 1 , 2 , TV 7 , respectively, along the diagonal. Since the matrix C consists of N 1 matrices 

*(3) 

K 7 ’ placed along the diagonal, its size is NpN^ x 5 NpN^. 
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Governing equations The governing equations for the interior and boundary subcells form a 
system of SNpN'y algebraic equations in the 5 NpN^ unknown coefficients This system of 

equations can be written symbolically as 

kT = t (59) 

where k is the structural thermal conductivity matrix. It contains information on the geometry 
and thermal conductivity of every subcell. The matrix k is made up of three submatrices as follows 


k = 


A 

B 

C 


(60) 


T is the thermal coefficient vector. It contains the unknown coefficients that describe the temper- 
ature field in all the subcells arranged as follows 


T = [T l5 ,T 7 , Tv 7 ] (61) 

where 

T 7 = [T (l7 \ , T^ 7) , T (jv ^ 7) ] (62) 

with T^ 7 ) defined previously in Eq. (52). Finally, t is the thermal force vector. It contains 
information on the boundary conditions specified either in terms of temperatures or heat fluxes. 

In summary, the elimination of generic cells in the user-friendly formulation has resulted in a 
system of equations which is easier to visualize and implement in a computer code. Most impor- 
tantly, the new formulation based on subcells as the basic building blocks sets the stage for an 
efficient reformulation presented in the next section. 

3 Reformulation 

The user-friendly formulation of HOTFGM-2D, which was presented in the preceding section, 
illustrates the feasibility of eliminating the distinction between generic cells and subcells. It should 
be noted that the simplified geometric model played an important role in the entire analytical 
process. In other words, the basic building block defined by the subcell itself without a generic 
cell, which was used to construct the heterogeneous composite, made the indices simpler and the 
derivation process more straightforward. More importantly, formulating the higher-order theory 
solely in terms of subcells, which play the role of generic subvolumes, reduces it to a standard 
thermal boundary-value problem on a discretized domain, thereby facilitating the search for a 
more efficient solution method. 

As discussed by Pindera (1991), the local/global stiffness matrix method is an efficient algorithm 
for solving mixed boundary-value problems in composite mechanics. This method is particularly 
well suited to problems involving layered composite materials or structures that require satisfaction 
of both continuity of tractions and displacements in mechanical problems, or temperature and 
heat flux in thermal problems, along common interfaces. The method thus far has been applied 
to layered media composed of continuous layers in Cartesian coordinates, or circular concentric 
rings/cylinders in polar coordinates. Examples of the method’s application include axisymmetric 
and axial shear problems involving multiple concentric cylinders, Pindera et al. (1993) and Williams 
and Pindera (1997), composite tube problems, Salzar et al. (1996), a layered composite cylinder 
subjected to diametral loading, Davison et al. (1994), and contact problems, Pindera and Lane 
(1993) and Urquhart and Pindera (1994). The local/global stiffness matrix method is based on 
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a systematic construction of a global stiffness matrix (mechanical problem) or global conductivity 
matrix (thermal problem) for the entire composite structure in terms of local stiffness matrices 
(mechanical problem) or local conductivity matrices (thermal problem) of individual layers. Herein, 
we extend this method to HOTFGM- 2 D where the volume discretization involves subvolumes of 
finite dimensions in the X2 — £3 plane defined by the subcells (/?y). Since we only discuss the 
thermal problem here, we will call it henceforth the local/global conductivity matrix method. 

In the local/global conductivity matrix method, the local conductivity matrix provides a con- 
nection between the surface- averaged heat flux and the corresponding surface- averaged temperature 
components at top, bottom, left and right surfaces of a given subcell. The assembly of local conduc- 
tivity matrices into the global conductivity matrix is based on the direct enforcement of continuity 
conditions along the subcell interfaces, and the interfacial temperatures are considered as the basic 
unknown variables in the reformulation. Therefore, for a volume discretized into NpN^y subcells, 
the total number of interfacial surface- averaged temperatures will be (Np + 1 )AT 7 + ( N . 7 + 1 )Np, 
which includes both the interior subcell interfaces and the external surfaces of boundary subcells. 

For a given (( 3 j) subcell, the local conductivity matrix can be expressed as follows, 


Q 2 

( fry ) 





( fry ) 

rr 1 — 

(£7) 

" K 11 

^12 

^13 

^14 " 

T 2 

Qi 


^21 

22 

^23 

^24 


1 2 


—Q 3 


^31 

^32 

33 

^34 




. Qi . 


_ ^41 

^42 

^43 

^44 


L ^3 J 



It relates the surface- averaged heat flux components — (Q^)^ 7 \ ~ (O3 )^ 7 \ ( O3" )^ 7 ^ on 

the left, right, top and bottom surfaces of a given subcell to their corresponding surface- averaged 
temperatures (T 2 - )^ 7 ), (T^)^ 7 ), ( T^)^\ Figure 8. The superscripts + and - refer to the 

corresponding heat flux or temperature acting on the right /top surface and left /bottom surface of 
a given subcell. The subscripts 2 and 3 indicate the X2 and x% directions, respectively. The minus 
signs in front of the heat flux components Q % , O3 indicate that the direction of the heat flux vector 
coincides with the direction of the unit normal vector associated with the left and bottom surfaces 
of the (/ 3 y) subcell. 



-fi 2 

f ~2 


Q + 3 ,t + 2 


subceU /3y 


e 2 + 

T + 2 


-Q3X2 


Figure 8. Definition of surface-averaged temperatures and heat fluxes. 
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Using the Fourier’s heat conduction law for orthotropic materials, Eqs. (33)- (34), 


% 


(07) 


= -k, 


(07) 


,97407) 


Ox 


( 0 ) ’ 


% 


(07) 


= -k. 


(07) 


Q 7 1 (07) 


Ox 


(7) 


the definitions of surface- averaged heat fluxes and temperatures in Eq. (63) are: 


and 


(Q+) ( 07) = 1 f ll/2 _fc(07) 
(7 7-^/2 


Q 7 1 (07) 


^( 0 ) i 4 /3) =^/2 


dx 


™(7) 


(Q-) ( 07) = 1 _ fc (07) 

<7 7—^/2 

ft, 3/2 


97^(07) 


<9^0) l 4 /3) =- ft /3/2 


dx. 


~(7) 


r ,2 -*r 8r,M 

hpJ-hpl 2 


~( 0 ) 


_( 7 ) l4 7 ) =E/2 


I r ^ 7 -?(^ 

I 4 ^=-Z t /2 a;E 2 


(^+)(07) = I f T (07) | 7^(7) 

2 ' U J-hl 2 =^/2 3 


t 7 7_z 7 /2 


1 /U/2 


(f-) ( 07) = I J 

™ m =iL 


Tl 


(07) 


^(7) 


7 J -I'll 2 
hp / 2 

h P J —hp/2 3 

ft/3/2 




=09) 


T (07) I 7 ~ 

Jq l 4 '» ) =i 7 /2 ax 2 


T, 


(07) 


U0) 


0 4-^/2 


l4 7j «U/2 ^2 


(64) 

(65) 

( 66 ) 

(67) 

( 68 ) 

(69) 

(70) 

(71) 


The detailed derivation of the elements of the local conductivity matrix in Eq. (63) will be outlined 
in the following section. 

Imposition of continuity of heat flux and temperature along common subcell interfaces, together 
with the external boundary conditions, gives rise to a system of equations in the unknown interfacial 
temperatures along the £2 and £3 directions, respectively. The continuity of interfacial heat flux is 
ensured by requiring that the sum of the heat fluxes acting on the interface separating (/?, 7 ) and 
(/? + 1 , 7 ) subcells and the interface separating (/?, 7) and (/?, 7 + 1 ) subcells be zero, i.e., 

( Q+) ( 0,7 ) + (_Q-)(0+ 1 ,7) = 0, (72) 

(Q+) ( 0,7) + ( _Q-)(0,7+ 1) = 0 (73) 


where (3 = 1 , 2 , ....Np — 1 and 7 = 1 , 2 , ....iV 7 — 1 . The thermal continuity is directly enforced by 
setting the surface- averaged temperatures along interfaces separating adjacent subcells to common 
values, 


(T 2 +)(0’ 7) 

(f+)(0,7) 


( T -)(0+ 1 , 7 ) _ j40+ 1 »7) 




(74) 

(75) 
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The system of equations is constructed by applying Eqs. (72)- (75) in conjunction with Eq. (63) 
to each interface, starting first with the left boundary in the first row of subcells where the boundary 
condition is prescribed and then progressing in the x 2 direction to each interface separating adjacent 
interior subcells, and then finally the right boundary with the prescribed boundary condition. This 
procedure is repeated along every row of subcells. We then repeat the same procedure in the X 3 
direction, starting at the bottom of the first column and proceeding up, and then repeating it 
for each subsequent column. At each interface, the thermal continuity expressed in terms of the 
common interfacial temperatures defined by Eqs. (74) and (75) is directly enforced in the heat 
flux continuity conditions given by Eqs. (72) and (73) through the use of Eq. (63). Actually, the 
assembly of the global conductivity matrix for the entire structure is carried out by superposing 
local conductivity matrices of individual subcells along X 2 and X 3 directions, respectively, in the 
manner that will be described in Section 3.3. 

Since the continuity of interfacial temperatures is enforced directly through Eqs. (74) and 
(75) along x 2 and X 3 directions, respectively, the above procedure results in the elimination of a 
substantial number of redundant interfacial continuity equations which are retained in the original 
formulation. Thus the reformulation simplifies the computational procedure for thermal boundary- 
value problems by reducing the size of the global conductivity matrix, which can be large for a large 
number of subcells. Furthermore, the global conductivity matrix can be generated systematically 
and automatically, thereby facilitating numerical implementation for any number of subcells. 

The above presentation illustrates the general formulation process of the local/global conduc- 
tivity matrix method in solving thermal boundary-value problems. In the following, the detailed 
derivation of how to construct the local conductivity matrix based on the generic cell-free formula- 
tion of HOTFGM-2D problem is outlined. It is should be noted that herein the thermal conductivity 
coefficients k and k 3^ are assumed to vary with local subcell coordinates. The same problem 
but with constant thermal conductivity coefficients was investigated by Bansal (2002). 


3.1 Local Conductivity Matrix: Variable Thermal Conductivity 

In the present work, the thermal conductivity coefficients k^^ and k 3^ in each (/3y) subcell are 
assumed to vary linearly with the local subcell coordinates x^ and x 3^ as follows 


k 

k 


(07) 

2 

(07) 

3 


- fr(07) , J07) A0) , h, (07 ) ^ ( 7 ) 

— ^20 a^ 22 rioQ ju q 


v 23 


— J07) , J07 ) f (P) , d,fa)Jn) 

regQ Ay q q Ju q 


v 32 


'33 


(76) 

(77) 


We note that k 23 ^ and k 32^ cannot be completely arbitrary as they do not explicitly appear in 
the surface- averaged heat fluxes and the volume- averaged heat conduction equation derived in the 
sequel. This is a direct result of the averaging technique employed in the higher-order theory. 
Therefore, we must impose the following constraint for consistency: k 23 ^ = k 32^ =0. In other 
words, we obtain one-dimensional variation of k^^ and k 3^ along the respective directions, if 
linear thermal conductivity variation is assumed. 

As before, the temperature field in the subcell (/?y) is approximated by a second-order expansion 
in the local coordinates x^\x^ in the following way (see Eq. (31)) 


7-407) 7-407) 1 A0) 7^(07) 1 ^(7)7-407) 1 _(0) 2 _ ^0 \ 7^(07) 1 .I^q^G) 2 

1 ^ X 2 1 (10) ' x 3 x (01) ^oW x 2 A ) 1 (20) oV 0X 3 


-( 00 ) 


( 10 ) 


-( 01 ) 


( 20 ) 


^7 \ 77(07) 
4 ; (02) 
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Using Eqs. ( 68 )-( 71 ), the surface- aver aged temperatures (T^)^ 7 ), (T 2 (T 3 are 

expressed in terms of the subcell temperature microvariables T^n ^ as follows 


( r r+\(P'y) 'T'iP'y) . 7 ) , P_rp((3 7 ) 

^2 J _i (00) 2 (10) + ^ ^ (20) 

( r v~\(/3^f) rpi.Pl) rri({3 7 ) | hp rp(/3y) 

v^2 J - i roo) 2 ( 10 ) + 71 r20) 


fi_ rri((3 7 ) 
. -*■ /on\ 


/7"H-\(/?7) rp(P*y) , ^7 rp((3~i) _|_ ^7 n^{Pl) 

' 3 I - ^ ( 00 ) + o UOl) + /( i f 02 ) 


( r r~\{Pl) rpifll) ^7 rfi{P / y) 

d 3 i _i (00) 2 ( 01 ) 


S T (W 
4 (02) 


Adding and subtracting the above first two equations involving (T^)^ 7 ) and (T 2 )^ 7 ^, we obtain 
expressions for and in terms of (Tr^)^ 7 ), (T 2 - )^ 7 ) and The matrix form of these 

two expressions become 


Similarly, and can be expressed in terms of (T s ( T ^)(^ 7 ) and T^j by applying 

the same technique to the above last two equations, 


The next step is to evaluate the surface-averaged heat flux components in terms of the mi- 
crovariables Integrating Eqs. ( 64 )-( 67 ) with k 2 ^ and k 3^ given by Eqs. ( 76 ) and ( 77 ), 

the surface-averaged heat fluxes (Q^)^ 7 \ (O2 (' Qt )^\ (O3 )^ 7 ^ are obtained in terms of the 

temperature microvariables T^\ as follows 

1 (mn) 


(f)+\{Pl) - rAfh)(_lSfh) _ n P h (fh)\ I rAfh) ( -^P utf 7) 

W2 J — 2 (10) l ^20 2 ft 22 ' J (20) ' 2 ^ 

(%-)<« = 

(n+\ (£7) — k ^ 

W3J — J (01) v %) 2^33 ' J (02) ' 2 ^30 

( q 3 -)( SiI = r <M ( _ fc <M + >, + - 


7 l(W\ I 71(^7)/ 6 h ljfil) 

2 Av 33 ; -f ( 02 ) ^ 2 ^30 

[lu(Pl)\ , t ,(/? 7 ) / 37 , (/ 3 7 ) 
K33 J T -f ( 02 ) l 2 K 30 


_ MpjJJhth 

a ^22 J 

Q/,2 

6n Ptm\ 


^ ^33 j 


7)) 

. /VQQ 


The above four equations can be rewritten in matrix form as follows 


-Q 2 1 < ‘’ 71 


^20 — koo ~r~ 


— k20 - ^22- 


Shpkzo 

2 

3 / 1 / 3^20 

2 


NASA/CR— 2002-21 1910 


27 



and 


' -<2 S _ ‘ 

(07) 

Qs J 



7 7 ^7 3Z 7 * 30 3Z 7 

*30 - *33 — I- — r *33 


u 


4 2 

7 7 ^7 3Z 7 * 3 o 3/ 7 

-* 30~*33 — — * 33 


(07) 


T ( 


( 01 ) 

r (02) 


(/ 3 'y) 


(81) 


2 2 

Substituting for t44 tV?) and T,4(l using Eqs. (78) and (79), the above equations become 


and 


Qs 

Qt 


(07) 


Q 2 

(07) 

o; 4fc 2 o , 2 fc 20 

2*22 , *22 7 

hp hp 

(07) 

' 4 ' 

(07) 

+ 

6^20 u 
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hp 
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7 2*20 oz 4*20 
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fys hp 


L 1 2 J 

6 * 2 ° xU 
+ 3*22 
hp 


1 (07) 
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-( 00 ) 


(82) 


2^33 - 
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&33 - 
-2/C33 


4^30 


4/^30 


f '7 J 


(07) 


tf 

T+ 

J 3 


(07) 


+ 


6fc 3 Q 

67730 

l 


- 3^33 

- 3^33 


7 


(07) 


T 


(07) 


( 00 ) 


(83) 


The last step is to develop a relationship between the surface-averaged temperatures and the 
volume-averaged temperature This relationship is obtained from the volume-averaged heat 

conduction equation 


f h/2 f h ^ /2 ( dq^_ + dq 3 

'-b/2 J-hp/2 y dxlf '* 


dx 


(07) \ 

(7) ) dx^dx^ = 0 


(84) 


where and q^ 1} are given by Eqs. (33) and (34). Since the thermal conductivity coefficients 


IPi) 


*2^ and k^ Jy) are functions of the coordinates x and x ^ y) , the above equation can be developed 
by using the chain rule as 

IJ2 h p / 2 d d 2 T {^) a 2 T(^), (/37 ). (j 9) (7) 

+ „_2 rm ^2 + „_iv> + „_2 iv> ^3 )d x 2 dx z — 0 (^) 


,(07) 


( 0 ) 


,( 7 ) 


Iry /2 h/3 /‘l 


dx)!p dx ^ <9^2 ^ 




Upon determining the derivatives of temperature and thermal conductivity, respectively, the above 
equation becomes 


I'y / 2 hpj 2 

f [ \b^T^ 4 - ( o*^ 7 ^ 4 - fr^ 7 )\ 4 - 

/ / L^22 1 (10) ^°^(20)V ZAi 22 x 2 ^^20 J ^ 


Substituting for T( 10 ) 
equation becomes 


—L y/2 —h/3 /2 

k(07) rp(Pj) | orn(Plf) 
o3 ^ (01) + ^(02) 


(/ A T (20) ^), and T( 


( 01 ) 


(2kf z l) x^ ) + k { 3 P)\dx^ ] dx^ = 0 (86) 

^ 7 ), T(o 2)^ 7 ^ using Eqs. (78) and (79), the above 
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k. 


(Pi) 

22 


[(T+)^) - (T-)^)] + [2(T 2 - )^ + 2(T 2 +)( /3 T 

J07) _ _ olW _ _ , a , 
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(00) J 


+ 


'7 


( 87 ) 


Therefore, the expression for T^ 7 ^ can be obtained from the above equation in terms of (T 2 )(^\ 
(f+)(Pl) 7 (Tg - )^), (T 3 + )^) as follows 
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( 88 ) 


Upon substituting the above expression for into the expressions for (Qjj - )^ 7 ), (Q 2 )^ 7 \ 

(Q 3 ~)^ 7 \ (Q 3 )^ 7 ^ given by Eqs. (82) and (83), the local conductivity matrix k^ 7 ) is obtained 
which directly relates the surface- averaged heat fluxes to the corresponding surface-averaged tem- 
peratures. These expressions can be written in matrix form as follows 
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(89) 


with the individual local conductivity matrix elements given below in terms of the subcell thermal 


conductivity coefficients and geometry, upon setting A = 12(/4 q 7 V^| + C’//?) for 
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3.2 Constant Thermal Conductivity 

When the coefficients associated with the local coordinates a: [7 and ai 3 7 ' ) in the expressions for 
and A'7 7) are set to zero, i.e., let A:^ 7 * = A-j/g 7 * = A'-^ 7 " 1 = A: 33 7 ' ) = 0 in Eqs. (76) and (77), constant 
thermal conductivities are obtained, i.e., = A:^ 7 * and A'7 7) = A: 3 q 7) ■ 

Replacing A^q 7 * with k^\ A' 3 q 7) with k^^and letting A'^ 77 = A' 33 7) = 0 in the local conductivity 
matrix k,^ 7 ^ with variable thermal conductivity in Eq. (89) , the new local conductivity matrix with 
constant thermal conductivities is obtained. 
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(07) _ „(07) _ 2 ^3 


(07) 


,(07) 


v 34 


k 43 — 


0 


-( 1 -^ 


2fe 


(07) ' 
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This is the same local conductivity matrix as that obtained by Bansal ( 2002 ), which confirms 
that the present development reduces correctly. 


3.3 Global Conductivity Matrix Assembly 

As mentioned above, the assembly of the global conductivity matrix is carried out by applying 
heat flux continuity conditions, Eq. (72) and (73), at each interface along x 2 and X 3 directions, 
respectively, together with the temperature continuity condition defined by Eqs. (74) and (75), 
starting at the boundary with prescribed boundary conditions. 


3.3.1 Temperature continuity conditions 

From the temperature continuity, we know that the surface-averaged temperature (T^)^ 7 ^ of the 
(/?, 7 ) subcell is the same as (T 2 _ )^ +1,7 ) of the (/3+1, 7 ) subcell defined in Eq. (74). The same holds 
true in the X 3 direction, Eq. (75). These conditions are applied f3 = 1, ..., Np — 1 and 7 = 1 ,..., N 1 — 1 
times at the subcell interfaces in the X 2 and X 3 directions, respectively, in the heat flux continuity 
conditions described below. Thus, altogether, the number of unknown surface- averaged interfacial 
temperatures in the interior is (Np — 1 )AT 7 + (N 7 — 1 )Np. 


3.3.2 Heat flux continuity conditions 

The continuity of interfacial surfaced-averaged heat flux along each interface requires that the sum 
of the heat fluxes across the interfaces separating adjacent subcells be zero, i.e., Eq. (72) in the X 2 
direction and Eq. (73) in the X 3 direction, respectively. Expressing these two equations in terms 
of the surface- averaged temperatures by using Eq. (63) yields for the interface between (/?, 7) and 
(/? + 1 , 7 ) subcells, 


4 / 4 7 ) f r r— 


V 21 


(Tr) m + k 42 ’^ ) ( 3 ' 2 + ) < ^ 1 + '4J 7 ) (?3"> ( '’ ,) + 




) (0+1,7) 


+47 1,7) (7 + ) (/3+1,7) + «i3 +1,7) (7“) (/3+1,7) + +1 ’ 7 ) (t 3 +) (/3+1 ’ 7) = 0 


(91) 


while for the interface between (/?, 7) and (/?, 7 + 1 ) subcells we get, 

47 ) (tr) (/37) + 42 ,7) (7 + ) (/?7) + 4J 7) (tr) (/?7) + «44 ,7) (7 + ) (/37) + 4i ,7+1) (7“) (/3,7+1) 
+4J 7+1) (7) (/3,7+1) + 4T +1) (7“) (/3,7+1) + 4J 7+1) 7 3 + ) (/3 ’ 7+1) = 0 ( 92 ) 

Employing the temperature continuity conditions, Eqs. (74) and (75), the above equations become 
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k 41 1 2 


, k .(A7)t=’(/ 3+1,7) , v/ ^, , i i-jrpxi^, I 1 -V | ... 

“I" ^42 ^2 “r ^31 ^2 “r ^32 

, /JAt) i (/3,7+l)\+>(/3,7+l) | (/3, 7+1) ^(/3, 7+2) _ ^ 
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43 ^3 


( 94 ) 


Thus, these two equations provide us with (Np — 1)TV 7 + (TV 7 — 1)7V^ equations for the unknown 
(Np — 1 )AT 7 + (iV 7 — l)Np common interfacial surface-averaged temperatures. 

3.3.3 Boundary conditions 

At external boundaries, the remaining 2(Np + N^) equations are obtained from the imposed bound- 
ary conditions specified either in terms of applied temperature or heat flux as given below 


7 ">(/M) I 
1 3 l® 3 =- 


or 


and 


or 


Qf l) 


T, 


OW 


\X3 = l N, 


_ _L [ hl! /2 T W dx W 
V2 _ ^+/2 “ 2 

t/2 2 
/2 ~h a J- h ./2 Ur 2 


q(P,N y) 


rhp/2 

l P J-h p /2 

rh f 3/2 
¥ 4-^/2 

dhp/2 

T’P J-hp/2 
1 _ Mi 0 w 

\X3=lN~y/2 L / Qtop aX 2 


(95) 

(96) 

(97) 

(98) 


+ 1-^/2 

at (/?, 1) and (/?, TV 7 ) subcells. Similar reasoning holds for subcells ( 1 , 7 ), and (TV#, 7) where the 
applied boundary conditions are given by 
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3.3.4 Governing equations 

The general form of the final system of equations can be written symbolically as 

Q = kT 

which can be decomposed for convenience as follows 
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. Q3 _ 
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where the heat flux vectors Q2 and Q3 are defined as 


Q 2 = 


— rnW 


Q 3 = [Q^,...,Qf,.. 


M^] 


( 105 ) 


where 


q( 7) = r nC 1 ^) 

Qf = [Q: 


n^+^h 


[Q 2 ’ IJ i 0, Q\ 

r\ o,...,o, Qf^«>] 


The vectors T2 and T3 contain the unknown surface- averaged temperatures at the subcell 
interfaces and the outer edges of the composite and are given by 


where 
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The global thermal conductivity matrix k is obtained by assembling the local conductivity 
matrices as explained above. It contains information on the geometry and thermal conductivities 
of the individual NpN^ subcells. The general structure of the global conductivity submatrices 
k 11, K12, k 2 i 5 k 22 is given below 
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The elements of the submatrices A^\ B^’ 7 \ and B^’ 7 ^ are given below 
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where 7 = 1,2, N . 7 . 



7 1) 4f } 0 0 0 0 0 0 
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3.4 Reformulated vs Unreformulated HOTFGM 

The motivation for the above outlined reformulation of HOTFGM is the reduction in the number of 
simultaneous equations whose solution determines the temperature field in the individual subcells. 
The reduction is a direct result of enforcing the thermal continuity at the common subcell interfaces 
directly. This is possible because the fundamental unknowns in the reformulation are the surface- 
averaged interfacial temperatures. In the original formulation, the size of the global conductivity 
matrix is 5 NpN 1 x 5 while in the reformulation version we have (2 NpN 1 + Np + N . 7 ) x 
(2 NpNry + Np + AT 7 ). Figure 9 compares the size of the global thermal conductivity matrices in 
the two versions. For discretized volumes with the same number of subcells in the two directions, 
Np = A^ 7 , the reduction can be as much as about 60%. 
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Figure 9. Reduction in the size of the global conductivity matrix due to reformulation. 
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4 Analytical Solution 


In this section, an analytical solution methodology is developed for determining temperature distri- 
butions in homogeneous and layered plates subjected to steady-state thermal loading. The solutions 
will be employed in validating both the user-friendly and reformulated versions of the higher-order 
theory. First, we develop a solution for a single strip or plate subjected to thermal boundary 
conditions, and then generalize it to layered plates consisting of an arbitrary number of strips. 

4.1 Single Plate/Strip Case 

Consider the problem of determining the temperature distribution in a rectangular plate whose 
length is a and height is b in the X 2 — £3 plane, Figure 10, subjected to steady-state thermal 
boundary conditions to be specified later. The thermal conductivity of the plate’s material is 
constant. Therefore, the governing equation for the temperature distribution within the plate in 
this problem is the Laplace equation 

V 2 T = 0 (107) 

We solve this equation by the separation-of- variables method, seeking a solution of the form 

T(x 2 , xz) = X (x 2 )Y ( 2 : 3 ) (108) 



I* ; ^ 


Figure 10. Problem definition. 


Substituting Eq. (108) into Eq. (107) and separating the variables gives 


so that 


and 


CN 

II 

1 

II 

* 1 * 


(109) 

x" -k 2 X = 0 


( 110 ) 

y" + k 2 Y = 0 


( 111 ) 

X(x2 ) ~ { A ^ x 2, 

1 C cosh kx 2 + D sinh kx 2 , 

k = 0 

k^0 

( 112 ) 
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Y(x 3 ) = 


E + Fx 3, 

1 G cos kx 3 + H sin kx 3, 


k = 0 
k ± 0 


(113) 


Because the Laplace equation V 2 T = 0 is linear, we can use superposition and combine the k = 0 
and k^O product solutions as follows 


T{x 2 , xs) = (A + Bx 2 )(E + TT 3 ) + ((7 cosh kx 2 + D sinh fcx 2 )(Gcos /CX 3 + H sinkx 3 ) (114) 

The coefficients A, , H depend on the particular boundary conditions applied at the plate’s edges. 

We consider two cases below: 

Case (i) 

The first case produces an internal temperature field that has no axes of symmetry when a/ 6 , 
and when a = b the temperature field is symmetric about the diagonal connecting the upper left 
and the lower right corners. The corresponding boundary conditions are 


T( 0 , 2 : 3 ) = 100 , T(a,x 3) = 200 

T(x 2 , 0 ) = 200 , T(x 2 , b) = 100 

Application of the boundary conditions at the bottom face £3 = 0, where T = 200, yields 


(. A + Bx 2 )E + (C cosh kx 2 + D sinh kx 2 )G = 200 


from which we obtain 

AE = 200, B = 0, G = 0 

Therefore Eq. (114) becomes 

T(# 2 , £ 3 ) = 200 + J ^3 + (C cosh/c £2 + T> sinh £^ 2 ) sin fcx 3 (115) 

where AF has been combined into J, CH into C\ and DH into D' for convenience. 

Next, application of the boundary condition at the top face x$ = 6 , where T = 100, yields 

200 + Jb + (C cosh kx 2 + D sinh kx 2 ) sin kb = 100 


from which we obtain 

J = —100/6, and sinfc 6 = 0 

Hence, 

k = mn/b (m — 1 , 2 , ) 

Therefore, Eq.(115) becomes 


T(x 2 ,x 3 ) = 200 


100^3 


+ ^(C m cosh m ™ 2 + sinh 


m7TX2 , . rmrxs 
1 sm ■ 


6 6 6 7 6 

n=l 

Application of the boundary condition at the left face X 2 = 0, where T = 100, yields 

m7TX 3 


200-±^2 + £ c „ 

m=l 


sm 


6 


= 100 


(116) 


(117) 


TT17YX 3 

Solving for the unknown coefficients C' m using the orthogonality properties of sin — - — , we obtain 


2 f b 

c - = iL { 


100;E 3 • nrnxs 

— 100 ) sm — - — ax 3 


(118) 
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or 


r __200 

1717T 


( 119 ) 


Finally, application of the boundary condition at the right face x 2 = a, where T = 200, yields 
100x 3 


200 


+ £< 

m=l 


200 , hitt a _ . , rrma x . ttityx^ 

cosh — h D m smh — — ) sin — — ± = 200 

71171 0 0 0 


(120) 


from which the solution for Dm can be obtained in the form 


Dm. = 


7 . _ , mna N , 
6smh( — - — ) 


7 6 T00x 3 . m 77X3 200 _ rriTva . 2 1 x 177 x 3 . 7 

/ ( — - — sm — 1 cosh — - — sin — - — )ax 3 (121) 

J 0 0 0 17177 00 


or 


Dm = 


r / rm7a. , 

200 [cosh ( — — ) — cosm7rJ 


ttity sinh( 


rmra. 


( 122 ) 


Thus, the solution to Eq. (107) is given by Eq. (116), with the coefficients C m and D m computed 
using Eqs. (119) and (122). 

Case (ii) 

If the boundary conditions are changed to produce temperature distribution which is symmetric 
about the X3 = 6/2 line, i.e., 

T(0,x 3 ) = 100, T(a, X3) = 200 
T(x 2 ,0) = 100 , T(x 2 , 6) = 100 

then we obtain a simpler form of solution. First, applying the boundary condition T(x 2 , 0) = 100, 
we get 

T(x 2, xs) = 100 + Jx 3 + (C cosh kx 2 + D sinh fcx 2 ) sin kx 3 (123) 

where AE = 100, B = 0, G = 0, and CH has been combined into C\ DH into D\ and AF into J. 
Next, application of the boundary condition T(x 2 , 6) = 100 in the above equation produces 


T(x 2 , 6) = 100 + Jx 3 + (C cosh kx 2 + D sinh fcx 2 ) sin kb = 100 
from which we obtain 

J = 0, sin kb = 0 

Therefore, the solution becomes 

00 

T(x 2, £3) = 100 + yyOn cosh sinh 


(124) 


n— 1 


rrmx2 x . rmix^ 
— ; — sm ■ 


6 


(125) 


After application of the boundary conditions T(0, X3) = 100 and T(a, x 3) = 200, the unknown 


uiTvy 

coefficients C m , are obtained, by using orthogonality properties of sin — - — , in the form 
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Dm = 


sinh ■ 


C m = 0 


inn • m7nr 3 , 

100 sm — 7 — ax 3 = 


400 


6 


m7r sinh 


rriTva 


(m = 1,3,5,...) 


(126) 

(127) 


6 6 

Thus, the solution in this case is given by Eqs. (125), (126), and (127). 
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4.2 N/r strip Case 

In the N^-strip problem, the rectangular plate, which has length a and height 6, is divided into N# 

strips along the x 2 direction. Therefore, every strip has height b and width ap, where (3 = 1, , Np, 

Figure 11. Since the boundary conditions on the top and bottom surface are given in terms of 
constant temperature, the analytical solution can be developed first from Eqs. (116) or (125) of 
the previous section depending on the applied boundary conditions. 



x 3 


(p) 


,'V 


(ap/2 ? h/2; 






/3 th strip 


Figure 11. N^-strip case problem definition. 

Case (i) 

Assume the boundary conditions on the top and bottom faces are T top = 200 and = 100. 
Then the temperature field in the arbitrary /3th strip is given by Eq. (116) as follows, 


T W(xf\ x 3 ) = 200 - + ]T (C<?> cosh + D W sillh !^3_) 


100^3 


mirx 


03) 


m nx 


03) 


b 


m— 1 


o . rriTYXs 
' sm — - — 
b 


However, it should be noted that this solution is obtained in the coordinate system with the 
origin located at the left-bottom corner. In order to obtain the solution in the coordinate system 
with the origin in the plate’s center, we use the following coordinate transformation, Figure 11. 

= x^ — a/ 3/2 => x^ = x^ + ap / 2 

xs = x^ — b/2 => xs = xs J rb/2 

Therefore, the solution of the Laplace’s equation in the translated coordinate system becomes 


= £ [Cg cosh + a * /2) + pW sinh + “‘ i/2) ] 

m=l 


gin m7r(x 3 + 6/2) + 2Q() _ 100(x 3 + 6/2) 
b b 


(128) 


Case (ii) 

If the boundary conditions are changed to produce temperature distribution which is symmetric 
about the x% = 6/2 line, i.e., 

Ttop = Tjjottom = 100 


NAS A/CR— 2002-2 11910 


39 


the temperature field in the arbitrary /3th strip is given by Eq. (125) as follows 

rOT(4»,x 3 ) = 100 + f; (COT cosh =£1 + D W slnh =£ ) sln = 

m=l 

Upon transforming the above solution from the coordinate system with the origin located at the 
left-bottom corner to that with the origin in the plate’s center, the solution of the Laplace’s equation 
with the symmetry boundary condition becomes, 

rW(4» * 3 ) = 100 + f) [COT cosh + + sinh TOT <y+'W2) 1 x 

m=l 

. m7r(x3 + 5/2) 

sin — — - 

b 

In both cases, the solution for the unknown coefficients Cm\ is obtained by constructing 
a system of equations which satisfies the heat flux and temperature continuity conditions at the 
interfaces separating the individual layers and the remaining boundary conditions. In order to 
make this as efficient as possible for large number of layers, we formulate the solution procedure 
using the local/global conductivity matrix approach. 

In the local/global conductivity matrix approach, the unknown coefficients cffl and D are 
the key elements in constructing the local conductivity matrix. They are expressed and calculated 
in terms of the harmonics of common interfacial temperatures between adjacent strips and subse- 
quently related to the corresponding interfacial heat flux harmonics through the local conductivity 
matrix. The detailed derivation is outlined in the following section. 



4.2.1 Construction of the local conductivity matrix 

In order to construct the local conductivity matrix, we start with the heat flux constitutive equation 
in the x 2 direction 

= ~ k t3-^py 0 = 1, (130) 

dx 2 

Substituting the expression for the temperature field given by Eq. (128) or (129), the above equation 
becomes 


,OT = 


m— 1 


COT sinh +-W 2 > + d OT cosh "”<4^+ <W?) 


. m7r(x 3 + 6/2) 

sin 


The determination of the coefficients cffl and D for case (i) and case (ii) is outlined below 


( 0 ) 


(131) 


Case (i) 

In order to express the coefficients cffl and D^} in terms of the surface temperatures on 
the vertical boundaries of the [3th strip, Eq. (128) is evaluated at (ap/2,x$) and (—ap/ 2,x 3 ), 
respectively, yielding 

2 , * 3 ) = 200 - 1W( Y h M + f; [COT cosh ^ + DOT 8 i„h =2] sin "*’ r( ^ + 6/2) 

m= 1 

(132) 
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or 


[T^(a p / 2,x 3 ) — 200 + 


100(S 3 + 6/2) 
6 


m 


jr [CW COS h =2 + D</> sinh = 2 ] S in 

m=l 

(133) 


and 


or 


TW(-a e /2,x 3 ) = 200 - 100 fe + 6 /2) + £ c fc» sin m*(S3 + 6/2) 

m=l 

[T<« (- 0 jJ /2, a*) - 200 + + = £ C W sin ” 1 ’ rtf3 i) + b/2) 

m— 1 


(134) 

(135) 


The corresponding interfacial temperature harmonics at the left (—ap/ 2, £ 3 ) and right (a^/2, £ 3 ) 
boundary for an arbitrary (3th strip can be obtained from Eqs. 133) and (135), respectively, using 
orthogonality relations 


= 


and 


6/ 2 
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— 6/2 

6/2 
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sin 02:3 


LXJ 

5 : [cff corf. ^2 + Dff sinh 


rrmap . 2 ^tt (^3 + V 2 ) 


sm 


dx 3 


m=l 


— 6/2 

6/2 

I 

-6/2 


(T 2 -)(f) = / T^\—ap/ 2 , X 3 ) — 200 + 


100(x 3 + 6/2) 


. mirixz + 6 / 2 ) 
sm ax 3 


f;cff>sm 3m ’ r <Y i,/2) ^3 


m=l 


Performing the required integrations, the above two equations become 

P?)£> = ^(Crcosh = + D M S mh=2: 

PPL® = |ctf> 

Upon writing the above two equations in matrix form, we obtain 


' T 2 ‘ 

09) 

T+ 


L z _ 

m 


6 


0 


6 rrmaR 6 rmrag 

- cosh — — * - - smh — — ! - 

2 6 2 6 


-Dm. 


l(/3) 


(136) 

(137) 


(138) 


Inverting the above system, the coefficients cffl and D are expressed in terms of the interfacial 
temperature harmonics as follows 
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D m 
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b rriirag b rrnrag 

- cosh — — ! - - sinh — — ! - 
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1 _1 

' t 2 - ■ 


L ^ 2 J 
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( 139 ) 


Similarly, evaluating the heat flux at the left {—ap/ 2, #3) and right (ap/2, £3) boundary for the /3th 
strip, respectively, yields 

5f«W2,*3) = = + 


4 ^ <-"V 2 - *3) = k/3 - r D m sin 


m=l 
00 


m7r . mir(x3 + 6/2) 


m=l 


b 


The interfacial heat flux harmonics for an arbitrary /3th strip are obtained from the above two 
equations by using orthogonality relations, 


(qW = 


6/2 

/ 

- 6/2 

6/2 


03 )/ /rj — \ mir(x3 + 6/2) 
(°/3 / 2 > x 3) sin cte 3 


(140) 


= / -fc„ £ (^Cff sinh ^ cosh =2) sin 2 ,rorfe + b/i) dx a 

m=l 


— (Q2 


- 6/2 

6/2 

/ 

- 6/2 

6/2 


4AN 0 /3/2, X3) sin 


. mir(x3 + 6/2) 


6 


dx 3 


/ ™ sin 2 ^(^3 + 6/2> dx, 

7/0 m=l 


(141) 


- 6/2 

The above two equations can be expressed in matrix form as follows 

0 mnk/3 

2 

mirka , miran mirka , rrmaR 
— ^ smh — — £■ — ^ cosh ■ p 


Q 2 

09) 

el . 

m 


6 


£>m 


03) 


(142) 


Case (ii) 


Similar to case (i), evaluating the surface temperatures on the vertical boundaries of the /3th 
strip i.e., Eq. (129) at (ap/ 2, £3) and (—ap/ 2,^3), respectively, yields 


T^\ap/ 2, x 3 ) = 100 + ^ [C^ cosh sinh 


m=l 


mirap , . m7r(x3 + 6/2) 
— sin — — - 


(143) 


or 


[TW („„/2, *3) - 100] m = f) [Off cosh ^ + gff sinh si, t/2) (144) 

m=l 
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and 


TW(-a,/2, 13) = 100 + f; Cff sin b/2) (145) 

m— 1 

or 

[tOT(-«,/ 2,S 3 ) - 100] m = f; C«/s (1 46) 

m=l 

The corresponding interfacial temperature harmonics at the left (—ap/ 2, S3) and right {ap/ 2, S3) 
boundary for an arbitrary /3th strip can be obtained from Eqs. (144) and (146), respectively, using 
orthogonality relations 
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Upon performing the above integrations, these two equations simplify to 

( T 2 + )m } = ^ (Cm ’ cosh } sinh ' 

(r 2 “)^ = ^ 

Writing the above two equations in matrix form, we obtain 
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6 rmrap 6 mirap 
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(147) 
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(149) 


which has the same form as Eq. (138) obtained for case (i) with the first kind of boundary 
conditions. Further, since the expression for the heat flux, Eq. (131), for both boundary conditions 
does not change, the relationship between the coefficients C and D and the interfacial heat 
flux harmonics, i.e., Eq. (142), is also the same. However, the specific values of the coefficients Cm\ 
D^} and the corresponding interfacial temperatures (T^)^, heat fluxes (Q^~)m\ (Qt)™) 

are not the same due to the different expressions for the solution of Laplace’s equation with different 
kinds of boundary conditions. 
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Nevertheless, in both cases, upon substituting the expression for the coefficients (7^ and D 
given by Eq. (139) into Eq. (142), carrying out the matrix multiplication and simplifying, we obtain 
the same relationship between the interfacial heat flux and interfacial temperature harmonics, i.e., 
the local conductivity matrix, 


where 

(^ll)rrP 
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(«2 !2)m^ 


Q 2 

03) 

K 11 

«12 

09) 

' t 2 ~ - 

Qt _ 

m 

_ «21 

« 22 

m 

L 1 2 J 


n 03) 


mirap 

mirks cosh — — — 
b_ 

. mirap 

b smh — - — 
b 

mirk/s 

_ mirap 

b smh — — — 
b 

mirks , miras 

— ^ sinh — — ^ + 

b b 

i u m7ra £ 

mirks cosh — — — 

6 

. mirap 
b smh — — — 
b 


2 mira/3 
kpmir cosh — - — 

I • u m7va (3 
b smh — - — — 
b 


(150) 


4.2.2 Global conductivity matrix assembly 

Imposition of the continuity of interfacial heat flux and interfacial temperature along common 
interfaces between adjacent strips, together with the external boundary conditions, gives rise to a 
system of equations in the unknown interfacial temperature harmonics. The continuity of interfacial 
heat flux is obtained by requiring the sum of the heat flux acting on the interface separating (3 and 
(3+1 strip to be zero, 

(Qt )^ } + (C 2 “)^ +1) = 0 P = 1, N P - 1 (151) 

whereas the continuity of interfacial temperature is given by 

(TW = (T 2 -)^ +1 > = if f3 = 1, ... .A^ - 1 (152) 

The system equations is constructed by using Eq. (151) at each interface along the X2 direc- 
tion, starting with the left boundary where the boundary condition is prescribed, together with 
application of the temperature continuity condition defined by Eq. (152), and proceeding to the 
right boundary. Following this procedure, we obtain the following system of equations for the mth 
harmonic of the interfacial temperatures 

left boundary: = — Q 1 ^ 

1st interface: ^ 21 T 1 ^ 1 + (^2 + ^ 11)^2 + ^ 12^2 = 0 

/ 3th interface: 2 + (ft^ 1 + ^n) T f~ l + = 0 

/ 3 + 1th interface: ^21 + (^22 + + 1^12 1 2 +1 = 0 

Nfs — 1 th interface: 1 T ^ f3 2 + + ^2 1 )T^ 1 + T(^ 9ht — 0 

right boundary: 1 + K^fT™ 9ht = Q™ 9ht 
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where Np is the number of strips, and the harmonic number m has been suppressed. The above 
equations can be expressed in matrix form for the rath harmonic as follows 
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The solution of the above system of equations for each interfacial temperature harmonic makes 
possible the determination of the coefficients C and D using Eq. (139). These, in turn, 
determine the temperature held in each strip of the plate using either Eq. (128) or (129). 


5 Mesh Sensitivity, Validation and Application Studies 

The user-friendly formulation of HOTFGM, developed in Section 2, dispenses with the two-level 
volume discretization employed in the original HOTFGM formulation wherein generic cells are fur- 
ther discretized into subcells. In the user-friendly formulation the graded microstructure is built 
up directly using subcells which play the role of fundamental subvolumes. This simplification in 
the volume discretization, which also simplified the derivation of the interfacial heat flux continuity 
conditions, provided the basis for the efficient reformulation outlined in Section 3. Therefore, the 
present section begins with the validation of the user-friendly formulation using the analytical solu- 
tions developed in Section 4. This is accomplished in two steps. First, we examine the convergence 
of the temperature field to the analytical solution as a function of the mesh discretization for the 
problem of a homogeneous plate subjected to the thermal boundary conditions discussed in Sec- 
tion 4. Once the convergence characteristics of the solution based on the user-friendly formulation 
are established, we then consider the problem of a piece-wise uniformly graded plate for which 
the analytical solution has also been derived in Section 4. This problem is subsequently re-solved 
using the reformulation of HOTFGM with variable thermal conductivity in order to demonstrate 
the reduction in the graded microstructure’s volume dicretization obtained by assuming linearly 
varying material properties. In the last example, a graded thermal barrier coating subjected to a 
through-thickness thermal gradient is considered, and the differences in the solutions based on the 
actual and homogenized microstructures are discussed. 

5.1 User-Friendly Formulation 

5.1.1 Mesh sensitivity study: homogeneous plate 

Let us consider a square plate, whose sides are unit length (i.e., a = b = 1), Figure 12, subjected 
to the following two sets of boundary conditions: 

(1) : Ti e f t = 100°G, T right = 200°G, T top = 100°G, T bottorn = 100° C 

(2) : Ti e f t = 100°G, T right = 200°C, T top = 100 o C, T boUo m = 200° C 
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4x4 subcells 



Figure 12. Problem definition and mesh discretization for a homogeneous plat( 
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The first set produces a temperature field that is symmetric about a horizontal line halfway across 
the plate’s height, whereas the second set produces a temperature field that is symmetric about 
a line connecting the upper left and lower right corners of the plate. Such temperature fields are 
sufficiently complicated to test the higher-order theory’s predictive capability. The plate’s thermal 
conductivity k has a constant value, i.e., k = 25 W/m — ° C. Since the plate is homogeneous, 
however, the temperature field does not depend on the conductivity under steady-state boundary 
conditions. This provided an additional check on the user-friendly formulation’s correctness as 
verified during preliminary numerical tests. 

The convergence behavior of the user-friendly formulation is investigated as a function of volume 
discretization by dividing the plate into 4 x 4, 12 x 12, 32 x 32 uniformly-spaced subcells, Figure 
12. The corresponding contour plots of the temperature distributions for the employed volume 
discretizations are compared with the analytical solution in Figure 13 for the first set of applied 
thermal boundary conditions. As is clearly seen in the contour plots, the temperature distribu- 
tions obtained from the user-friendly formulation are symmetric with respect to the horizontal line 
halfway across the plate’s height irrespective of the volume discretization, and gradually approach 
the analytical solution with increasing number of subcells. In the case of the 4x4 subcell mesh, 
we obtain only an approximate temperature field which, nevertheless, exhibits the basic character- 
istics of the analytical solution. Increasing the volume discretization to 12 x 12 subcells generates 
a temperature field which appears quite close to the analytical solution except at the upper and 
lower right corners where the boundary temperature distributions change suddenly, producing dis- 
continuities at these corners. Increasing the mesh to 32 x 32 subcells eliminates these local corner 
disturbances, producing smooth curve patterns in the temperature field which are visually the same 
as those obtained from the analytical solution. 

A more quantitive picture of the user-friendly formulation’s convergence behavior with mesh 
refinement is obtained by plotting the temperature distributions along specific cross sections. Three 
horizontal cross sections and three vertical cross sections have been chosen as illustrated by the 
arrows in Figure 12. The horizontal cross sections are situated at the elevations x% = 0.0, 0.25, 0.5, 
and the vertical cross sections are situated at X 2 = 0.5, 0.75, 1.0. The cross sections include the lower 
horizontal boundary defined by x$ = 0.0 and the right vertical boundary defined by X 2 = 1.0. The 
interior cross sections have been chosen such that they run along subcell interfaces. For example, 
the cross-section situated at xs = 0.25 runs along the horizontal interfaces that separate the sets of 
subcells (/? 1) and (/32) in the case of the 4x4 subcell mesh. This same cross section runs along the 
horizontal interfaces that separate the sets of subcells (/? 3) and (/34) for the 12 x 12 subcell mesh. 

As we know from the user-friendly formulation, the temperature distributions along the differ- 
ent interior cross sections can be obtained by using microvariables associated with the subcell on 
either side of the common interface. Furthermore, recall that the thermal and heat flux continuity 
conditions were applied in a surface-average sense. Therefore, the temperature distributions along 
a given cross section separating two sets of subcells, obtained using two different sets of microvari- 
ables associated with either set of subcells, will generally be different. However, these distributions 
should approach the actual temperature distributions with increasing mesh refinement. Further- 
more, the surface-averaged interfacial temperatures, which exhibit piece-wise uniform distributions, 
should approach the actual temperature distributions in a step-wise manner with increasing mesh 
refinement. 

In order to demonstrate the above behavior, three types of temperature distributions along 
the chosen subcell interfaces have been generated for the interior cross sections. In the following 
horizontal cross-section plots, we define the common interfacial surface-averaged temperature as T, 
the interfacial temperature obtained using the microvariables of the subcells below the interfaces 
as T + , and the interfacial temperature obtained sing the microvariables of the subcells above the 
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Figure 13. Temperature fields for the three mesh discretizations and first set of boundary 
conditions. Comparison with the analytical solution. 
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interfaces as T~ . Similar notation is used for the vertical cross-section plots. In these cases, T + 
indicates the interfacial temperature obtained using the microvariables of the subcells lying to 
the left of the subcell interfaces and T~ denotes the interfacial temperature obtained using the 
microvariables of the subcells lying to the right of the subcell interfaces. In the case of the boundary 
subcells, only two curves are given since there is only one set of subcells forming the boundary. 
As before, the designation T is used to denote the surface-averaged temperature distribution along 
both the horizontal and vertical boundaries. The second curve is obtained using the properties of 
the boundary subcells: for the lower horizontal boundary it is denoted by T~ and for the right 
vertical boundary it is denoted by T + for consistency with the established notation. 

The temperature distributions along the interior horizontal cross sections at the two elevations 
£3 = 0.25 and X 3 = 0.5 for the first set of boundary conditions are shown in Figure 14. For the 
4x4 subcell mesh, the temperature distributions calculated on the different sides of the subcell 
interfaces in the X 3 = 0.25 cross section shown in the left column of Figure 14 are quite different, 
with discontinuities or jumps evident as the subcell vertical interfaces are traversed along the X 2 
direction. The differences, however, gradually disappear with increasing mesh refinement, including 
the above-mentioned jumps. In the case of the 12 x 12 subcell mesh the differences are already quite 
small and barely visually discernible, while in the case of the 32 x 32 subcell mesh the differences 
cannot be identified visually. These distributions are graphically indistinguishable from the actual 
distributions obtained from the analytical solution (not shown). The surface- averaged temperature 
distributions also gradually converge to the actual temperature distributions in a step-wise manner. 

The corresponding distributions in the x% = 0.5 cross section are shown in the right column 
of Figure 14. The same observations with regard to the solution’s convergence with increasing 
mesh refinement hold in this cross section as in the preceding one, with one important difference 
however. Interestingly, we only see two curves in this cross section, namely the surface-averaged 
temperature distribution T and the two T + and T~ distributions superposed on each other. This 
is due to the symmetric boundary conditions employed in this example, which in turn produce a 
temperature field that is symmetric with respect to the X 3 = 0.5 cross section. This, in turn, ensures 
that the subcell microvariables on either side of the line of symmetry are the same, producing the 
same interfacial temperature distributions. This behavior is clearly captured by the user-friendly 
formulation, providing further evidence of its correctness. 

Figure 15 shows the temperature distributions along the vertical cross sections X 2 = 0.5 and 
X 2 = 0.75 as a function of the volume discretization for the first set of boundary conditions. In 
this case, the temperature distributions are not symmetric about any of the two cross sections, 
and consequently three distinct sets of interfacial temperature distributions, T, T + and T~ are 
present. Symmetry exists only about the x% = 0.5 cross section, and this is reflected in each of the 
individual distributions for each cross section (since the temperature distributions are given as a 
function of the x% coordinate) and for each volume discretization. As in the case of the horizontal 
cross sections, the difference between the T + and T~ distributions vanish with increasing mesh 
refinement (including the jumps at the subcell horizontal interfaces), while the T distribution 
approaches the actual distribution in a step-wise manner. 

The temperature distributions at the two external boundaries X 3 = 0 and X 2 = 1 are given in 
Figure 16 for the three volume discretizations and the first set of boundary conditions. From these 
figures, we conclude that the T~ temperature distributions for the lower horizontal boundary, shown 
in the left column of Figure 16, and the T + temperature distributions for the right vertical boundary, 
shown in the right column of Figure 16, gradually approach the applied boundary conditions. 
However, even with 32 x 32 subcell mesh, they exhibit some fluctuations at the lower and upper right 
corners of the plate. That is because the applied temperature distributions possess discontinuities 
at these corners, requiring further mesh discretization in the vicinity of these points. 
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Figure 14. Convergence of temperature distributions along the cross-sections X 3 = 0.25, 0.5. 
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Figure 15. Convergences of the temperature distributions along the cross-sections x 2 = 0.5,0.75. 
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Figure 16. Convergence of the temperature distributions along the edges x% = 0 and X 2 = 1. 
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The corresponding temperature contour plots and cross section distributions for the second set 
of boundary conditions are given in Figures 17 through 20. The full field contour distributions 
shown in Figure 17 illustrate temperature field symmetry with respect to the cross section passing 
through the upper left and lower right corners of the plate for all three volume discretizations. 
As in the preceding case, the results obtained from the user-friendly formulation of the higher- 
order theory converge to the actual temperature distribution obtained from the analytical solution 
with increasing mesh refinement. The full field temperature distribution obtained from the 32 x 32 
subcell mesh is visually indistinguishable from the distribution generated by the analytical solution. 
Similar observations with regard to the convergence behavior with increasing mesh refinement hold 
for the cross section temperature distributions along interfaces separating the interior subcells, 
presented in Figures 18 and 19, and along external boundaries presented in Figure 20. In this case, 
the temperature distributions are not symmetric about the chosen cross sections as expected, and 
therefore three distinct sets of interfacial distributions are evident in the interior, which converge 
with increasing mesh refinement at the same rate as before. 

In summary, as observed from the above figures, whatever interfaces are considered, for either 
set of the boundary conditions (Figures 13-15 or Figures 17-20), the differences between the tem- 
perature distributions based on the microvariables belonging to adjacent subcells decrease with 
increasing number of subcells, while the surface- averaged temperature distributions approach the 
actual distributions in a step-wise manner, with the step increments becoming increasingly smaller. 
In other words, the greater the mesh refinement, the better the point- wise convergence or accuracy 
to the actual distribution. 

5.1.2 Validation: discretely graded plate 

In this section, the results from the user-friendly formulation for a plate with discretely graded 
thermal conductivity in a piecewise uniform manner are compared with the corresponding analytical 
solution outlined in Section 4 for the N^-strip plate. The effect of thermal conductivity variation 
on the temperature field is also investigated. 

We consider the same boundary value problem discussed in the preceding section but limit the 
investigation to the first type of boundary conditions, i.e., the plate is subjected to a temperature of 
200°C on the right boundary and a temperature of 100°(7 on the remaining boundaries, Figure 21. 
The plate is divided into vertical strips and the thermal conductivity k (assuming isotropic thermal 
properties) is assigned a different constant value in each strip, producing a piecewise uniformly 
graded plate. The overall variation of the thermal conductivity is given by exponential functions 
in the x 2 coordinate. Two different thermal conductivity variations are employed. One is an 
exponentially increasing function of X 2 

k = 0.5e 5x2 

while the other is an exponentially decreasing function of X 2 

k = 60e~ 5x2 

These continuous functions are employed to produce piecewise uniform variations as follows. 
The plate is divided into 25 vertical strips of equal width and the above functions are used to 
calculate the thermal conductivity at the center of each strip, which is then assumed to remain 
constant in that strip, as shown in the two lower graphs of Figure 21. The analytical solutions 
for this problem are then obtained based on the N^-strip solution outlined in Section 4 using 
a 200 harmonic representation of the temperature field in each of the 25 strips. The correspond- 
ing user-friendly formulation solution is obtained by dividing the plate into 25 x 25 subcells of equal 
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Figure 17. Temperature fields for the three mesh discretizations and second set of boundary 
conditions. Comparison with the analytical solution. 
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Figure 18. Convergence of the temperature distributions along the cross-sections X 2 = 0.5,0. 
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Figure 19. Convergence of the temperature distribution along the cross-sections xs = 0.25, 0. 
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Figure 20. Convergence of the temperature distributions along the edges x% = 0 and ^2 = 1. 
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Figure 21. Problem definition and mesh discretization for a layered plate with piece- wise 

uniformly graded thermal conductivity. 
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Figure 22. Temperature field comparison between HOTFGM-2D and analytical solutions for the 
layered plate with exponentially increasing (top) and decreasing (bottom) piece-wise uniformly 

graded thermal conductivity. 
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dimensions, shown in the upper right schematic of Figure 21, and assigning the same values of 
thermal conductivity in each vertical column of subcells as in the case of the analytical solution. 
The above boundary conditions and manner of grading will produce a temperature field which is 
symmetric about the horizontal cross section halfway across the plate’s height. 

The comparison of the full field temperature distributions obtained from the analytical and 
higher-order theory solutions is given in Figure 22 for the exponentially increasing (top) and de- 
creasing (bottom) thermal conductivity variations. In both cases, excellent agreement is evident 
throughout almost the entire plate with the exception of isolated regions in the vicinity of the 
upper and lower right corners. These regions exhibit large temperature gradients caused by the 
temperature discontinuities at the corners themselves. Local mesh refinement is therefore needed 
to produce better correlation in these regions in the case of HOTFGM. 

As is also observed in Figure 22, when the thermal conductivity increases in the x 2 direction 
according to the function k = 0.5e 5x2 , the temperature field is shifted towards the left vertical 
boundary with respect to the temperature field in the homogeneous plate subjected to the same 
boundary conditions shown in Figure 13. Alternatively, when the thermal conductivity decreases 
in the X 2 direction according to the function k = GOe -53 ’ 2 , the temperature field is shifted towards 
the right vertical boundary. This can be explained by the Fourier’s heat conduction law, 

q\ Pl) = -k^djT^ 

where the temperature field is seen to be modulated by the spatial variation in thermal conductivity. 

5.2 Efficient Reformulation: Variable Thermal Conductivity 

In this section, the efficient reformulation of the higher-order theory described in Section 3 is em- 
ployed to generate the temperature field in a continuously graded plate whose thermal conductivity 
variations are described by the exponential functions of the preceding section. To mimic the ex- 
ponentially varying thermal conductivity in the X 2 direction, three volumetric discretizations were 
employed. In all cases the plate was discretized into 25 subcells along the x% direction, while along 
the X 2 direction the plate was discretized into 10 and 12 uniformly spaced, and 12 nonuniformly 
spaced subcells, resulting in two uniformly spaced 10 x 25 and 12 x 25 microstructures and one 
nonuniformly spaced 12 x 25 microstructure. The thermal conductivity k^ 1 ^ in each vertical col- 
umn of (/? 7 ) subcells was varied linearly with the local subcell coordinate x^ such that continuity 
was preserved from one column of subcells to another along the global X 2 coordinate. The values 
of the thermal conductivity at the subcell boundaries were calculated using the exponential func- 
tions employed in the preceding section, and subsequently employed in the construction of linear 
approximation of the thermal conductivity k^ 1 ^ in each subcell. Figure 23 illustrates the piece- 
wise linearly varying k^ 1 ^ in the manner that approximates the exponentially increasing thermal 
conductivity for the three volume discretizations. The thermal conductivity in the X 3 direction, 
kf l] , was kept constant in each column of subcells, but varied from one column of subcells to the 
next by taking the average value of the thermal conductivities at the subcell boundaries calculated 
from the exponential functions. The general expressions for variable thermal conductivities in the 
X 2 and X 3 directions, k^ 1 ^ and k^\ respectively, thus reduce to the following expressions, 

dPi) _ l.(P 7 ) 1 7 (/ 3 7 )-(/ 3 ) 

^2 — ^20 “ r " ^22 x 2 

iSM) — ju (/* y ) 

^3 — ^30 
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Figure 23. Thermal conductivity approximation by piece-wise linear variation. 

The full field temperature distributions obtained for the 10 x 25 subcell mesh with the exponen- 
tially increasing (top) and decreasing (bottom) thermal conductivities, generated using the efficient 
reformulation with locally linearly varying conductivities, are presented in the right column of Fig- 
ure 24. The corresponding predictions generated using the user-friendly formulation with 25 x 25 
subcell mesh of the preceding section are included in the left column of Figure 24 for comparison. 
The agreement is very good everywhere except in the upper and lower right corner regions where 
the temperature jumps occur. Increasing the mesh size to uniformly spaced 12 x 25 subcells im- 
proves the agreement in these corner regions, as seen in the left column of Figure 25. Retaining 
the same mesh size but preferentially increasing the number of subcells in the corner regions where 
large thermal gradients occur produces essentially the same temperature distributions in the entire 
linearly graded plate as those in the discretely graded plate, as seen in the right column of Figure 
25. This indicates that a smaller number of subcells with local linearly varying thermal conduc- 
tivity is required to capture the essential character of the temperature field in a discretely graded 
plate with piecewise uniform thermal conductivity. 
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Figure 24. Temperature field comparison between user-friendly formulation and efficient 

reformulation of HOTFGM-2D. 
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Figure 25. Temperature field comparison between evenly and nonevenly distributed thermal 
conductivity predicted by the efficient reformulation of HOTFGM-2D with linearly varying 

thermal conductivity. 
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5.3 A Thermal Barrier Coating Application 

In this section, the temperature field in a graded thermal barrier coating subjected to a through- 
thickness temperature gradient is investigated using the efficient reformulation of the higher-order 
theory. The zirconia coating is graded with metallic inclusions such that the volume fraction of the 
inclusion phase increases with increasing distance from the hot surface of the coating exposed to the 
temperature of 200° ( 7 , Figure 26. The cold surface where the metallic inclusion content is greatest 
is exposed to 100 °C. The coating’s graded microstructure, which is periodic along the vertical 
direction, is simulated by the zero heat flux boundary conditions on the top and bottom surfaces 
of the coating. The thermal conductivities of the zirconia coating and the metallic inclusions were 
taken to be 0.5 and 100 W/m — ° ( 7 , respectively. 




Figure 26. Actual microstructure and volume fraction distribution of the inclusion phase. 

The thermal analysis is conducted using the actual graded microstructure seen in Figure 26 
which has been discretized into 10 x 100 subcell mesh, and two additional discretely graded and 
linearly graded microstructures with local homogenization. The locally homogenized microstruc- 
tures were obtained from the actual microstructure by first determining the local metallic inclusion 
content in each of the nine regions into which the actual microstructure was divided in the manner 
shown in Figure 26. Subsequently, each region with its local inclusion content was modeled using 
HOTFGM as a strip containing 16 uniformly spaced metallic inclusions subjected to a known macro- 
scopic temperature gradient along its length and zero heat flux on the top and bottom surfaces, 
Figure 27. 
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Figure 27. Homogenized thermal conductivity vs fiber volume fraction. 


The effective thermal conductivity for such strip was calculated from the ratio of the average heat 
flux obtained from the solution of the specified boundary-value problem and the negative of the 
imposed average or macroscopic temperature gradient. The choice of 16 uniformly spaced metallic 
inclusions was dictated by the asymptotic behavior of the average thermal conductivity as a func- 
tion of the inclusion number. The results of the calculations are included in Figure 27 in the form 
of a graph of the macroscopic thermal conductivity as a function of the metallic inclusion volume 
fraction which corresponds to each of the regions shown in Figure 26. The resulting discretely and 
linearly graded in a piecewise manner thermal conductivities along the actual coating’s length are 
given in Figure 28. The piecewise linearly graded microstructure was obtained from the piecewise 
discretely graded microstructure by using the center points of adjacent subcells to define new sub- 
cells with linearly varying conductivities. Due to the homogenization of the microstructure and the 
applied boundary conditions which produce a one-dimensional temperature field, the homogenized 
microstructures were represented by 1 x 9 subcell meshes. 

The full field temperature distributions in the actual and homogenized graded coatings are 
given in Figure 29. The temperature distribution in the actual coating is clearly two-dimensional, 
with the disturbances created by the inclusion presence more pronounced near the upper and lower 
boundaries than in the interior. In contrast, the temperature fields in the homogenized coatings 
are clearly one-dimensional, with a slight shift present between the discretely and linearly graded 
homogenized microstructures. A more quantitative picture of the differences in the temperature dis- 
tributions among the three microstructures is obtained by comparing the cross section plots. These 
are presented in the top portion of Figure 30 with three cross section temperature distributions 
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Figure 28. Homogenized thermal conductivity as a function of position. 


in the actual coating and one in each of the two homogenized coatings. The staircase tempera- 
ture pattern observed in the actual coating is due to the pronounced differences in the thermal 
conductivities of the metallic inclusion and ceramic matrix phases. In contrast, the temperature 
distributions in the homogenized coatings are smoothly and continuously increasing functions of 
the coordinate along the coating’s length, with the differences between the two distributions visible 
near the hot surface. The bottom portion of Figure 30 compares the temperature distribution in 
the actual coating averaged across the coating’s width with the distributions in the homogenized 
coatings. In both cases, the temperature distributions in the actual coating, averaged or unaveraged 
across the coatings width, are above the temperature distribution in the discretely graded coating. 
This is not true in the case of the linearly graded coating wherein the temperature distribution is 
actually higher in some regions near the hot surface than the temperature distribution in the actual 
coating. 

Despite the fact that the temperature distributions in the two homogenized coatings are almost 
identical with the exception of the region near the hot surface, the average heat flux calculated 
in the homogenized coating with piecewise linear variation of the thermal conductivity is closer to 
the average heat flux in the actual coating than the averag flux in the discretely graded coating. 
These values can be used to calculate the macroscopic thermal conductivity of the coating across 
its thickness as shown below, 
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Figure 29. Temperature distributions in TBC with actual and homogenized microstructures. 
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Figure 30. Cross-section temperature distributions comparison. 
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where AT = 100° C is the temperature difference between the top and bottom surfaces of the 
coating and L = 100 is the coating’s thickness. Consequently, the effective thermal conductivity 
of the entire coating based on the linearly graded microstructure is closer to the effective thermal 
conductivity of the actual coating than the discretely homogenized coating. Thus it appears that 
the linear thermal conductivity variation is more appropriate for estimating the effective thermal 
conductivities of actual graded microstructures than the discrete thermal conductivity variation 
when the coating’s microstructure is relatively coarse. 


6 Conclusions and Future Work 

The goal of this research was to develop a generic cell-free, subcell-based formulation of the higher- 
order theory for functionally graded materials developed by Aboudi et al. (1999), which would 
then provide a basis for a more efficient reformulation using the local/global conductivity matrix 
approach. In order to extend the higher-order theory’s range of applicability, variable thermal 
conductivity capability at the subcell level was also incorporated into the efficient reformulation. 
This investigation was undertaken because the original formulation of HOTFGM is computationally 
intensive. This, in turn, limits the size of problems that can be analyzed due to the large number 
of equations required to mimic realistic microstructures of functionally graded materials. 

In order to validate both the user-friendly formulation and the reformulation of HOTFGM with 
variable thermal conductivity, analytical solutions to the steady-state heat equation were devel- 
oped for homogeneous and layered plates. The homogeneous plate solution was used to study the 
convergence behavior of the user-friendly formulation as a function of mesh discretization. The 
layered-plate solution was employed to investigate temperature distributions in discretely graded 
plates. Comparison between the results obtained from the user-friendly formulation and analytical 
solutions validates the simplified set of equations of the higher-order theory based on the more 
straightforward volume discretization which does not differentiate between generic cells and sub- 
cells. In particular, identical results were obtained for a homogeneous plate subjected to thermal 
boundary conditions that produced internal temperature fields with different symmetries. The so- 
lution’s convergence was demonstrated to be quite rapid. A piecewise uniformly graded plate with 
exponentially varying thermal conductivity was also analyzed and the results obtained from the 
higher-order theory coincided with the analytical solution for a layered plate. 

The user-friendly formulation is easy to understand and implement. More importantly, it pro- 
vides a basis for the development of an efficient reformulation based on the local/global conductivity 
matrix approach, which was also successfully accomplished, taking the locally variable thermal con- 
ductivity at the subcell level into account. The efficient reformulation has significantly decreased 
the number of equations. The reduction of the number of equations by as much as 60% for some 
problems has enhanced the theory’s capability, enabling solution of more practical and complicated 
problems. Comparison of temperature fields in piecewise uniformly graded and linearly graded 
cases, both of which approximate the same spatially varying thermal conductivity at the global 
scale, revealed that the same temperature distribution can be obtained with a smaller volume dis- 
cretization when the thermal conductivity is allowed to vary at the subcell level. Therefore, the 
efficient reformulation of the higher-order theory with variable thermal conductivity made possible 
the solution of a wide range of problems requiring large number of subcells. Thus, after transi- 
tion from the user-friendly formulation to reformulation, the temperature distribution prediction 
for arbitrary microstructures under many circumstances can be accurately accomplished in a more 
efficient manner. 
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However, this research was limited to steady-state heat conduction of graded materials in two 
dimensions. In order to further extend the theory’s range of applicability to practical problems, 
further research is required. The next step is to extend the presented efficient reformulation to 
enable the analysis of mechanical problems. As the first step, the two-dimensional elastic refor- 
mulation of the higher-order theory with linearly varying material properties within each subcell 
will be developed and then extended to three dimensions. This will complement the efficient ther- 
momechanical reformulation of the higher-order theory with constant subcell material properties 
recently completed by Bansal (2002). Further, in order to make full use of material properties and 
optimize structural design, it is necessary to incorporate viscoelastic, viscoplastic or plastic effects 
of material response into the theoretical framework. Therefore, viscoelastic, viscoplastic and plastic 
capabilities should be included in the reformulated higher-order theory. Fracture analysis is also an 
important area of FGMs. Therefore, the capability to analyze fracture mechanics problems should 
also be included in the reformulated higher-order theory. 
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